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ABSTRACT 


This  thesis  describes  the  study  and  development  of  a 
workable  inverse  scattering  method  for  imaging  and  identifi¬ 
cation  of  radar  targets.  The  space -time  integral  approach 
is  used  for  iterative  target  shape  reconstruction. 

Following  an  overview  of  transient  electromagnetics,  the 
integral  equation  is  applied  for  thin-wire  transient  response 
computation. 

The  analytical  time  domain  integral  equation  is  derived 
and  solved  numerically  for  general  conducting  bodies  of 
revolution.  Finally,  the  algorithm  for  an  inverse  scattering 
computer  solution  is  derived  and  tested  under  simulation  of 
physical  environments. 
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INTRODUCTION 


I . 


A.  OVERVIEW 

The  classification  of  radar  targets,  has  become  increasing¬ 
ly  important  in  the  past  two  decades.  In  military  applications, 
the  need  to  distinguish  between  friendly,  enemy,  and  neutral 
targets,  is  vital. 

In  most  radar  applications,  the  only  properties  that  are 
measured  are  range,  angle,  velocity  and  elevation.  The  methods 
that  are  used  by  radars  to  provide  some  classification  of  the 
target  generally  involve  the  examination  of  the  echo  signal. 

Some  of  these  techniques  can  be  summarized  briefly  as  follows: 

1.  High  range  resolution  radar,  which  uses  a  very  short 
pulse  that  can  provide  sufficient  range  resolution  to  obtain 
a  rough  profile  of  the  illuminated  part  of  the  target  shape. 

2.  Cross  section  fluctuations  analysis  which  can  provide 
information  on  the  size  of  the  target. 

3.  Synthetic  aperture  radar  that  gives  a  high  angular 
resolution.  This  method  is  typically  employed  in  ground 
mapping . 

4.  Analysis  of  the  backscattered  polarization  from  a  target 
can  provide  a  means  for  discrimination  between  types  of 
targets . 

5.  Use  of  nonlinear  effects  due  to  junction  of  metals. 

The  backscattered  spectrum  is  analyzed  and  checked  for 
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harmonics  appearing  in  the  received  signal,  that  are  typical 
for  certain  type  of  targets. 

6.  Doppler  modulation  by  propellers,  the  rotary  wing  in 
helicopters  or  even  hull  vibrations  in  armored  vehicles. 

None  of  the  methods  cited  can  be  used  to  obtain  all  the 
information  on  the  target  shape  and  size,  they  only  can  be 
used  to  discriminate  to  a  certain  extent  between  different 
types  of  targets.  Other  methods  currently  used  involve 
passive  reception  of  radiation  by  targets. 

Elint  receivers  or  IR  detectors  can  be  used  to  monitor 
radiation  from  targets,  and,  based  on  previous  intelligence, 
or  signature  catalog,  associate  a  type  of  radiation  to  a  type 
of  target. 

The  inverse  scattering  technique , that  is  the  subject  of 
study  in  this  thesis,  uses  the  basic  properties  of  the  transi¬ 
ent  response  to  an  impulse  excitation  (impulse  response)  in 
order  to  obtain  information  on  the  size,  shape  and  also  material 
composition  of  the  target. 

In  principle,  the  smoothed  impulse  response  can  be  ob¬ 
tained  by  measuring  the  backscattered  fields  at  all  frequen¬ 
cies  up  to  some  upper  limit.  In  most  cases  this  v/ould  be 
impractical  so,  instead,  the  target  is  illuminated  by  a  very 
short  pulse  (with  a  very  wide  band  frequency  spectrum  ranging 
from  the  low  frequencies  in  the  Rayleigh  region  up  to  higher 
resonance  region  frequencies) . 
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The  impulse  response  contains  all  the  information  needed 
to  characterize  the  target.  The  higher  frequencies  give  the 


information  about  the  structure.  This  is  related  to  the 
optical  theory  of  diffraction  which  is  a  good  approximation, 
under  certain  conditions,  for  shape  prediction  up  to  the 
shadow  boundary.  Physcial  optics,  however,  fails  to  predict 
the  phenomena  appearing  in  the  shadow  region,  such  as  the 
creeping  wave  which  is  due  to  lower  frequency  waves  travelling 
around  the  body . 

A  different  approach  that  has  been  used  to  retrieve  tar¬ 
get  configuration  information  is  related  to  the  singularity 
expansion  method  (S.E..M.).  This  method  seeks  to  classify 
scatterers  through  the  complex  residues  and  poles  of  their  time 
domain  signatures.  Briefly,  in  this  method  the  transient  E.M. 
response  of  a  conducting  body  is  characterized  as  a  series  of 
complex  exponentials  (or  damped  sinusoids):  [1],  [2] 

N  St 


where  R  is  the  amplitude  (residue)  of  each  mode  of  complex 
m 

frequency  s^.  The  true  function  has  a  complex  frequency 
representation  of  the  form: 


F(s) 


N 

I 


R 


ra 


1  s  -  s 
m=l  m 


The  poles  depend  only  on  the  geometry  (independently 
of  the  orientation) ,  and  the  residue  depends  on  both 
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excitation  and  geometry  [3].  Thus^once  the  poles  and  residues 
are  known,  they  provide  a  means  for  target  classification. 

The  problem  is  to  find  the  poles  and  residues. 

A  method  is  suggested  using  Prony's  algorithm  in  12)  for 
extraction  of  poles  and  residues  from  time-domain  signatures. 
The  method  has  been  tried  successfully  on  wire  geometries  and 
simple  shapes. 

Another  important  application  of  the  impulse  response  of 
a  target  is  the  study  of  its  radar  cross  section  (RCS) ,  since 
the  Fourier  transform  of  the  impulse  response,  (the  frequency 
response)  is  directly  related  to  the  RCS  by  the  following 
relationship : 

c  =  -  c^  iF(u)) 

where 

F((^)  =  /  F^(t)  dt 

—  oo 

and 

F^(t)  is  the  impulse  response  of  the  target. 

Finally  the  time-domain  analysis  of  the  impulse  response  of 
a  body  is  of  great  importcuice  in  the  study  of  the  electro¬ 
magnetic  pulse  (E.M.P.)  response  of  targets. 
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B.  PRACTICAL  CONSIDERATIONS 


The  practical  study  and  measurement  of  the  impulse  re¬ 
sponse  involves  the  building  of  a  time-domain  scattering  range. 
The  availability  of  generators  having  very  short  pulse-widths, 
fast  rise  time,  and  high  voltage  outputs,  together  with  fast 
sampling  oscilloscopes  and  small  computers  makes  the  task 
feasible . 

For  the  solution  of  the  inverse  scattering  problem  the 
time  domain  scattering  range  developed  at  the  Naval  Post¬ 
graduate  School  [4],  was  employed  to  measure  the  impulse 
response  and  demonstrate  the  validity  of  the  theories  developed 
in  this  thesis.  The  transient  response  data  measured  by  using 
pulse  generators  with  very  fast  rise  times  can  validate  the 
calculation  of  the  time  domain  integral  equation  solution. 

More  important,  it  provides  a  vehicle  for  testing  various  tar¬ 
get  imagery  and  I.D.  schemes,  as  well  as  investigating  broad¬ 
band  RCS  reduction. 

C.  SURVEY  OF  THE  LITERATURE 

One  of  the  first  efforts  to  determine  physical  properties 
of  electromagnetic  scatterers  from  time-domain  analysis  was 
due  to  Kennaugh  and  Cosgriff  [5] ,  who  employed  the  physical 
optics  approximation  to  calculate  the  approximate  backscattered 
impulse  response  of  a  flat  spheroid. 

A  summary  of  this  work  was  reported  by  Kennaugh  and  Moffatt 
[6] ,  who  extended  the  physical  optics  approximation  to  transient 
response  calculations.  At  this  stage  the  computation  of  the 
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impulse  response  was  approximated  by  using  the  physical  optics 
approach  or  the  inverse  transform  of  the  frequency  domain  solu¬ 
tion.  Then,  in  1968  ,  Bennett  derived  in  his  Ph.D.  the'.is  [71 
an  approach  for  solving  the  time  domain  integral  equation 
directly  by  using  a  time -stepping  method.  Bennett  used  a 
form  of  the  equation  called  the  magnetic  field  integral  equation 
(or  M.F.I.E.)  for  two  and  three-dimensional  surfaces. 

Sayre  and  Harrington  [8]  used  the  same  time-stepping  tech¬ 
nique  to  derive  the  transient  response  of  a  wire,  but  they 
employed  another  type  of  equation  called  the  electric  field 
integral  equation  (EFIE) . 

Additional  work  in  time-domain  electromagnetics  has  been 
reported  in  several  more  recent  papers  as,  for  example,  the 
report  by  Miller,  Poggio  and  Burke  [9]  on  time-domain  analysis 
of  thin-wires.  Summaries  and  review  of  the  time-domain  electro¬ 
magnetics  are  due  to  Mittra  (2]  and  [10]  and  Miller  [111. 

A  comprehensive  review  on  the  subject  and  its  application 
has  been  published  by  Bennett  [12],  which  gives  the  reader 
an  introduction  and  initiation  to  time-domain  electromagnetics. 
Literature  dealing  with  the  inverse  scattering  problem  is  due 
to  Young  [13]  and  Shubert,  Young  and  Moffatt  [14],  which  des¬ 
cribes  the  image  generation  of  a  target  from  the  ramp  response 
waveform. 

Moffatt  and  Mains  [15]  describe  a  method  of  detection  and 
discrimination  of  radar  targets  using  the  idea  of  ramp  response 
but  with  multiple  harmonic  radar  frequencies.  Finally ,  a  very 
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recent  report  by  Bennett  [16J  gives  the  general  approach  to 
the  time  domain  inverse  scattering  solution.  The  method  sug¬ 
gested  in  that  report  is  being  used  in  this  thesis  to  approach 
the  problem  of  transient  target  imaging.  This  list  of  efforts 
is  by  no  means  exhaustive,  as  many  other  authors  have  published 
papers  dealing  with  transient  electromagnetics  and  target 
imaging . 

D.  THESIS  OBJECTIVES 

The  objective  of  this  effort  is  to  develop  the  analytical 
and  numerical  solution  of  time-domain  inverse  scattering 
problem.  It  follows  the  development  of  the  time-domain 
scattering  range  at  the  NFS  [4]. 

The  main  effort  has  been  conducted  toward  the  formulation 
and  numerical  solution  of  the  time-domain  integral  equations 
for  conducting  bodies  of  revolution. 

A  computer  progreim  is  developed  to  solve  for  direct  and 
inverse  scattering.  The  algorithm  was  used  for  experimental 
display  of  the  shape  of  the  target  under  test. 
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II.  GENERAL  SOLUTION  OF  THE  TIME 
DOMAIN  INTEGRAL  EQUATIONS 

This  chapter  is  a  review  of  the  development  of  time  domain 
integral  equations,  their  solutions  and  the  associated  numeri¬ 
cal  considerations.  Detailed  development  of  the  equations 
can  be  found  in  [2,10,11]. 

A.  THE  INTEGRAL  EQUATIONS  (E.F.T.E.  and  M.F.I.E.) 

The  basis  for  derivation  of  the  equations  are  the  general 
time  dependent  forms  of  Maxwell's  equations: 

7  ^E(r,t)  =  -  ^  H(r,t) 

J  t  o 

7  xH(r,t)  =  ^  E{r,t)  +  J(r,t) 

ot  O 

7  •  E{r,t)  =  o(r,t)/c^ 

7  •  H^(r,t)  =  0 

Starting  from  these  equations,  and  using  the  expressions  for 
the  magnetic  vector  potential,  A,  and  the  electric  scalar 
potential,  4- »  which  are  related  to  the  fields  by 

H{r,t)  =  -i-  V  xA{r,t) 

E(r,t)  =  -  V  (p(r,t)  -  |Y(r,t) 
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one  can  obtain  the  expression  of  the  magnetic  field  integral 
equation  (M.F.I.E.)  for  a  perfect  conductor  [2] 


Jg(r,t)  =  2n  <H^^^(r,t) 


2'T  'gCOT  RS 


(2.1) 


and  the  electric  field  integral  equation  (EFIR) 


n<E-"(?,t)  = 


+  -  |:r)RldS' 
.2  R  Cot 


£  R 
o 


(2.2) 


In  equations  (2.1)  and  (2.2), 


n  is  the  unit  vector  normal  to  the  conductor  at 
the  observation  point  r. 

J  (r',r)  is  the  current  density  at  the  source  point  r' 
®  at  a  retarded  time  x  . 


:(r','')  is  the  charge  density  at  source  point  at  the 
retarded  time  i  . 


Here  :  is  given  by 


X  =  t  -  R/c 


and 


R  is  the  distance  from  the  observation  point  r  to 
the  source  point  i.e., 

R  =  r  -  r' 
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E  (r,t)  and  H  (r,t)  are  the  incident  electric  and  magnetic 
fields  at  the  observation  point.  Because  of  numerical  con¬ 
straints,  the  equation  most  suitable  for  closed  surfaces  is 
the  MFIE  equation  (2,1)  while  for  thin  wires  and  open  struc¬ 
tures  the  EFIE  equation  (2.2)  is  used. 

Equation  (2.1)  states  that  the  current  induced  on  a  point 

—  -•inc 

r  of  a  body  at  a  time  t  is  due  to  the  incident  field  H  at 

that  point  and  the  sum  of  scattered  fields  there  from  earlier, 

more  distant  locations  of  the  body.  The  interactions  between 

the  samples  of  the  body  are  displaced  in  time  by  an  amount 

equal  to  the  transit  time  of  the  fields  between  them. 

In  the  integral,  the  observation  point  (R  =  0)  is  excluued 
since  this  "self-term"  is  included  in  the  incident  field  por¬ 
tion  of  the  equation.  This  determines  the  general  solution 
for  the  MFIE  equation  which  is  solved  directly  by  marching 
on  in  time.  The  solution  is  started  at  a  certain  time  (t  =  0) , 
and  due  to  causality  the  integral  part  of  the  equation  will 
be  zero  since  there  are  no  past  currents.  Then,  as  t  increases, 
the  integral  is  solved  with  previous  currents  already  computed. 
This  technique  is  a  time-stepping  procedure  for  solving  initial 
value  problems. 

For  the  EFIE,  equation  (2.2),  the  method  is  slightly  differ¬ 
ent  but  uses  the  same  time  stepping  procedure.  The  solution 
of  this  equation  is  applied  to  the  problem  of  transient 
scattering  by  a  thin  wire,  in  the  next  chapter. 

The  solution  of  equations  (2.1)  and  (2.2)  gives  the  cur¬ 
rents  induced  on  the  body  by  the  incident  field.  Once  the 
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currents  are  known,  it  is  a  straightforward  step  to  compute 
the  scattered  fields.  For  the  far-field  case  the  field  is 
given  by: 


,t) 

o 


4TTr  c 
o  S 


^  J{r* ,t: 

0  T 


“is 


(2.31 


where  a.^  is  the  unit  vector  from  the  source  r'  to  the  far 
observation  point  and  r^  is  the  separation  distance.  J^tr',!) 
is  the  current  at  the  source  point  at  a  previous  (retarded) 
time  T  =  t  -  r  /c. 


B.  NUMERICAL  SOLUTION  AND  CONSIDERATIONS 

The  solution  of  the  integral  equations  implies  finding 
the  currents  induced  on  the  body  once  illuminated  by  an  inci¬ 
dent  field.  Since  the  solution  is  to  be  carried  out  numerically, 
the  first  step  is  to  represent  the  currents  by  their  space  and 

time  samples  J^j(s,t)  in  a  discrete  manner.  J^j(s,t)  is  the 

t  h  t  h 

current  at  the  i^  segment  LS  (or  patch)  at  the  j  time 

increment  AT.  The  total  current  will  be  approximated  in  terms 

of  all  samples  at  all  times  and  spatial  segments  via  the  basis 

function  expansion: 


J  (r ,  t) 


I  I  J  (s,t)  U(t  )  V(s  ) 
i=l  j=l  ^  ^ 


(2.4) 


where  Ng  is  the  number  of  space  samples  on  the  body  while  N^ 
is  the  total  number  of  time  increments,  and  the  pulse  basis 
functions  are: 
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U(tj)  = 


V(si) 


1  at  the  center  of  the  jth  time  interval 


0  elsewhere 


1  at  the  center  of  the  ith  space  patch 


0  elsewhere. 


The  space  time  sampling  is  shown  in  figure  2.1,  The 
sampling  of  the  currents  already  suggests  that  the  body  is 
divided  into  segments  (or  patches)  AS  and  the  time  is  sampled 
at  AT  intervals,  while  the  total  observation  time  goes  up  to 
'<  AT. 

It  is  very  important  at  this  point,  before  attempting  to 
solve  this  equation,  to  define  the  basic  criteria  for  proper 
selection  of  the  space  and  time  sampling  as  well  as  their 
relationship  to  the  shape  of  the  exciting  pulse  and  the  body 


The  incident  field  is  chosen  here  to  be  represented  by  a 
gaussian  impulse  due  to  its  rapidly  decreasing  features  and 
also  because  most  of  the  short  pulse  generators  will  produce 
such  a  waveform. 

The  gaussian  impulse  is  given  by 


g(t)  =  e 


-A^t^ 


Its  amplitude  falls  to  1/10  of  its  maximum  value  at 
tmax  “  1.5/A.  Its  frequency  spectrum  is  given  by  [2] 

G(f)  =  F{g(t)  }  =  e  ^  . 
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utaticn 


2] 


This  response  reduces  to  1/10  of  its  maximum  value  at 

f  1  0.5A.  The  time  sampling  AT  must  exceed  the  Shannon 
ma  X  r 

rate  to  avoid  aliasing  problems  if  FFT's  are  employed  later, 
i  .  e .  , 


AT 

min 


1 


2f 


max 


1 

A 


Usually  AT  will  be  taken  smaller. 


AT 


1 


5f 


max 


1 

2.5A 


(2.5) 


The  spatial  sampling  AS  must  be  chosen  such  that 

AS  ^  AT  xc  (2.6) 

otherwise  there  will  be  coupling  between  the  currents  of 
adjacent  patches.  Violation  of  (2.6)  will  also  mean  that  the 
stepping  procedure  for  solving  the  equations,  which  assumes 
that  during  one  time  interval  the  currents  at  one  patch  are 
not  influenced  by  adjacent  currents,  will  no  longer  be  satis¬ 
fied,  thus  producing  an  incorrect  solution. 

Another  consideration  is  that  the  space  sample  points  must 
be  close  enough  to  resolve  adequately  the  spatial  variation 
of  the  incident  field  over  the  body.  From  Eq.  (2.5)  and  (2.6) 
we  have: 


AS 


-  5f 


max 
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This  latter  relationship  imposes  some  minimum  ratio  of  the 
size  L  of  the  body  to  the  minimimi  wavelength  (or  equivalently 
the  minimum  pulse  width) .  The  ratio  L/A  determines  the 
resolution  and  the  usefulness  of  the  time  domain  analysis, 
and  it  must  be  kept  high  enough  to  maintain  adequate  resolution. 
If  the  body  is  to  be  represented  with  good  fidelity  with  a 
certain  minimum  number  of  samples  we  will  need 


L 


'  min 
5 


This  means  that  if  A  .  increases  (i.e.,  wider  pulse)  L  must 

min  f  ' 

increase  or  equivalently  for  a  small  size  target  X  must 

ram 

decrease  (i.e.,  shorter  pulse). 

Having  presented  the  equations,  their  general  solution  and 
the  requirements  for  the  various  parameters  involved  we  can 
move  on  to  the  implementation. 

C.  APPLICATION  TO  THE  TIME-DOMAIN  SCATTERING  BY  A  THIN-WIRE 
The  first  application  that  was  studied  was  the  computation 
of  the  impulse  response  of  a  straight  thin  wire.  This  example 
was  chosen  since  it  gives  a  good  understanding  of  the  physical 
processes  and  the  numerical  techniques  used  in  the  solution 
of  time  domain  integral  equations.  The  basis  for  this  compu¬ 
tation  is  given  by  Sayre  and  Harrington  in  [8] . 

The  problem  is  formulated  as  follows:  a  thin-wire  of 
length  L  and  radius  a  aligned  with  the  z  axis  is  illuminated 
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In  these  equations  I{z',t)  and  q(z',T)  are  the  current  and 
charge  at  source  point  z'  at  time  t  =  t  -  R/c ,  while  R  is 
the  distance  between  the  source  and  observation  points.  The 
geometry  is  described  in  figure  2.2. 

Using  the  numerical  formulation  expressed  in  [8],  and 
applying  the  time  stepping  procedure,  a  computer  program  was 
written  (Appendix  D)  to  compute  the  scattered  field  from  a 
thin  wire  of  length  1  meter  and  radius  0.8  mm.  The  wire  was 
divided  into  40  segments  and  the  sampling  time  AT  =  AZ/c  =  L/40c. 
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The  incident  field  was  taken  to  be  similar  to  the  pulse 


transmitted  by  the 
range,  with  linear 
incidence  normally 


generator  of  the  NPS  transient  scattering 
polarization  parallel  to  the  wire  and 
on  it. 


E 


inc 


(t-t, 


max 


2 


o  •  10® 


-L 

sec 


t  =  18*1T  (about  half  tra.isit  time  along 

m3X  .  , 

the  wire) 


The  program  computes  the  currents  at  each  spatial  segment  with 
up  to  400  time  samples.  This  corresponds  to  10  transit  times 
across  the  wire.  The  current  at  the  center  of  the  wire  is 
given  in  figure  2.3.  The  backscattered  field  is  given  in 
figure  2.4. 

No  significant  change  was  obtained  for  20  segments  and 
200  time  increments.  Comparison  of  this  computation  with  lab 
measurement  was  done  by  Hammond  [4]  and  very  good  agreement 
was  found. 

Similarly,  comparison  of  the  computed  results  found  in 

Q 

[11]  with  a  ratio  of  L/a  =  0.00667  and  A  =  3.25-10  gave 
exact  agreement.  The  backscattered  field  for  this  case  is 
shown  in  figure  2.5. 
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Q  [11]  Figure  58 


Figure  2.5.  Backscattered  field  from  a  thin-wire 
(A  =  3.25.1o9  and  a/L  =  0.0067) 


III.  TIME  DOMAIN  RESPONSE  OF  PERFECTLY 
CONDUCTING  BODY  OF  REVOLUTION 

A.  INTRODUCTION 

This  section  considers  the  solution  of  the  time  domain 
response  for  a  body  of  revolution,  when  the  exciting  wave  is 
propagated  along  the  axis  of  symmetry  of  the  body. 

The  use  of  the  axial  symmetry  of  the  body  will  result  in 
considerable  simplification  of  the  computation  of  both  the 
induced  currents  on  the  body  surface  and  the  field  scattered 
from  the  body.  The  algorithm  developed  in  this  section  will 
be  used  later  for  the  solution  of  the  inverse  scattering 
problem. 

B.  THE  MAGNETIC  FIELD  INTEGRAL  EQUATION 

In  [7]  Bennett  derives  the  integro-dif ferential  equation 
for  three  dimensional  scatterers: 


J  (r ,  t) 


2a  X  H  (r , t) 
n 


where : 


*  Jb  )la„  + 


2'  a''^" 


RC  9t' 


(3.1) 


I  =  t  -  R/C  (retarded  time) 
r  is  the  observation  point 
r'  is  the  source  (or  integration)  point 


30 


R  =  :  r  -  r ' 

A  is  the  surface  of  the  scatterer 

H^(r,t)  is  the  incident  magnetic  field  at  the  observer's 
point. 


The  geometry  for  the  body  of  revolution  is  given  in  figure 

3 .1. 


For  the  actual  problem  the  incident  magnetic  field  is 
propagating  along  the  z-axis  (axis  of  rotation  of  the  body) , 
and  is  directed  towards  the  positive  x-axis.  The  body  lies 
completely  to  the  right  of  the  positive  z-axis,  and  is  defined 
by  the  continuous  contour  function,  p(z) ,  which  is  the  radius 
of  the  body  at  a  point  z. 

The  prime  coordinates  denote  the  source  or  integration 
point,  while  the  unprimed  denote  the  observation  point  for 
which  the  current  is  desired.  At  each  point  on  the  body,  there 


is  defined  the  triad  of  orthogonal  unit  vectors, 
where ; 


a  :  unit  vector  normal  to  the  surface 
n 

a  ;  unit  vector  tangent  to  the  surface  and  directed 
^  towards  the  positive  <})  direction 

a^ ;  unit  vector  tangent  to  the  surface  and  directed 
along  the  curve  representing  the  longitudinal 
contour  of  the  body. 


For  the  following  derivation  we  will  define  the  angle  a  (or 

a')  as  that  sustained  by  a  (or  a')  with  the  z  axis,  oi  (or  a') 

s  s 

being  positive  when  a^  (or  a^)  points  outward  from  the  z  axis 
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(or  j ' ) 


and  negative  if  it  points  towards  it.  The  angle  p 
is  the  angle  defined  by  (or  .  ’)  with  the  (x,y)  plane. 
Equation  (3.1)  states  that  the  current  density  at  the  obser¬ 
vation  point  r  at  time  t  is  due  to: 

1.  The  current  induced  by  the  incident  field  at  (r,t). 

This  part  represents  the  physical  optics  approximation, 

2a  xH^{r,t). 
n 

2.  The  current  due  to  the  flow  of  currents  at  all  other 
points  of  the  body  at  a  retardated  time  t  =  t  -  R/C,  where 
R  is  the  distance  from  the  integration  point  to  the  oberva- 
tion  point. 

According  to  the  geometry  described  in  figure  3 . 1  we  can  write 

R  =  +p'^  -  2.;p’  cos  +  (z’ 

The  current  on  the  body  will  flow  only  in  the  s  and  $  direc¬ 
tions  and  can  be  expressed  (for  both  primed  and  unprimed 
coordinates) ,  as 

J(r,t)  =  J„  (p  ,  z,  t,(p) -a^  +  (p  ,z,t,()))  .a, 
s  S  <J)  i 

where  J  (p,z,t,4i)  and  J  (p,z,t,(j))  are  the  scalar  values  of 

S  p 

the  current  components  and  are  functions  of  the  position  on 
the  scatterer.  We  will  denote  them  simply  by  J  and  J.  so 

S  'p 

that 
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Eq.  (3.1)  can  be  written  in  terms  of  its  components  in  the 
s  and  3  direction,  i.e., 


Jg (r , t) 


[2a^  ^HNr,t)].  +  F  ^  ^A' 


2-  ■  C  3t'  ^“n 


(3.2) 


J^(r,t)  =  [2a^  <  (r,t)  ].  i  |^]  [a„  ^  j  ■  x  r].  dA' 


2-n  C  3t^  ‘  n 

A  K 


(3.3) 


In  order  to  obtain  the  scalar  expressions  of  the  two  last 
equations  we  have  to  perform  the  vector  multiplications: 

This  is  done  in  Appendix  A  and  as  a  result,  equations  (3.2) 
and  (3.3)  become: 

A  R 

(3.4) 

J  =  -2H^(z,t)cos  >  sin  a+.5^  ^  +  i]  [ J ’  U  sin  6 +J IF, ] dA’ 

,  L.  <7  T  K  S  9  ^ 

(3.5) 


where : 


=  I  z  '  -  z)  cos  i3  sin  J  +  (p  -  '  cos  3  )  cos  a  ' 

V  =  z '  -  z 

F^  =  (  P  cos  3  -  o' )  cos  a  +  sin  ot  cos  3  (  z  '  -  z) 

'J  =  ( z  '  -  z )  sin  a  •  sin  a'  +p  cos  a  sin  a '  -  p  '  sin  a  •  cos  a  ' 
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We  note  that  the  incident  tie  Id  is  dependent  only  on  z  and 
t  (independnet  of  i  due  to  s'_,Tnme try )  .  We  also  note  that: 

dA’  =  .  •  di’  ds’  =  •  dp  ds’ 

The  integrals  in  Eq.  (3.4)  and  (3.5)  can  then  be  divided  into 
two  parts:  integration  over  ti  from  0  to  2n,  and  integration 
over  ds'  from  0  to  S  (S  being  the  contour  curve  as  shown  in 
figure  3.2). 


The  next  step  is  to  obtain  further  simplication  of  the  inte¬ 
gral  equations  (3.4)  and  (3.5). 

The  currents  in  these  equations  are  functions  of  their 
position  on  the  body: 

~  (  z,p  ,(&  ,  t) 

But  we  remark  that  the  currents  due  to  the  incident  field  have 
a  sinusoidal  variation  with  ,  as  follows: 

J,  =  H^(z,t)  cos  (1)  •  sin  a 

'P 
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Jg  =  ( z , t)  sin  i 


so  that  we  can  assume  that  the  currents  will  have  the  same 
respective  relationships  with  $ . 

J  (z,p,t,i)  =  J  (z,^,t)sin  $  =  J  -sin  4)  (3.6) 

s  s  s 

J  (z,t,t,4))  =  ( z ,  p  ,  t)  cos  i  =  J  -cos  o  (3.7) 

This  assumption  is  demonstrated  in  Appendix  B,  using  Maxwell’s 
equation . 

After  substitution  of  equations  (3.6)  and  (3.7)  into  the 
integrals  (3.4)  and  (3.5)  ,  the  integrals  over  the  angle  g  can 
be  simplified  as  shown  in  Appendix  C,  so  that  the  scalar 
integral  equations  reduce  to  their  final  expressions: 


[J'F,  cos  3  +Vsin^6J']d6 
si  9 


(3.8) 


S  27t 

■2H^  (z,t)  sin  a  +  ^  /  p'ds'  / 


[J^  U  sin^S  +  cos  g  ^ 


(3.9) 


in  which  the  currents  J  and  (J*  and  J')  are  only  functions 

S  vf  S  4) 

of  p,  z,  and  t  (p',z',r)  and  independent  of  the  angle  i . 

The  significance  here  is  that  the  surface  distribution  of 
the  unknown  currents  is  reduced  to  a  one-dimensional  variation 
of  current  along  the  generating  arc  of  the  surface  of  revolution. 
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This  will  allow  considerable  reduction  in  computer  time  and 
storage  when  carrying  out  the  numerical  solution  of  the 
equations . 

Once  the  currents  along  the  contour  of  the  body  are  known, 
the  currents  at  any  other  point  (at  some  time  t  and  some  z 
position)  will  be  given  by 

Jg(p,z,t,j)  =  Jg  ( >: ,  z,  t)  •  sin  5 

J^(p,z,t,f)  =  J  (p,z,t)-cos  4). 

C.  FAR -ZONE  SCATTERED  FIELD 

In  most  cases  we  are  interested  in  the  far  scattered 
field.  At  a  large  distance  from  the  body  the  scattered  field 
is  given  by: 


H®(r,t)  =  /:^-[J(r',T)  xa^^JdA'  (3.8) 

o  A  ' 

where 


:  is  the  retardated  time  t  -  R'/C. 

R'  is  the  path  difference  between  the  integration 
point  to  the  observer  and  the  origin  to  the 
observer. 

J(r',T)  is  the  current  at  the  surface  integration  point 
retarded  in  time,  while  a„  is  the  unit  vector 
pointing  from  the  source  point  towards  the  obser¬ 
vation  point.  The  geometry  is  described  in 
figure  3.3. 

r  is  the  distance  from  the  origin  to  the  observation 
°  point. 
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A  is  the  projection  of  the  source  point  on  the 
plane  (x,z) 


Figure  3.3.  Geometry  for  the  far-zone 
scattered  field 
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In  the  following,  the  scattered  field  will  be  derived 
for  an  arbitrary  elevation  angle  .  Due  to  symmetry  with 
r  there  is  no  loss  of  generality  if  we  observe  the  scattered 
field  in  the  plane  (x,z)  where  ;  =  0. 

According  to  the  geometry: 

R'  =  -[z'  cos  '  COSO'  sin  6] 

T  =  t  -  Ff/C. 

R'  is  negative  whenever  the  observation  point  is  closer  to  the 
source  than  the  origin,  and  positive  otherwise. 

a_  =  sin  a,  +  cos  a 
R  X  z 

Performing  the  vector  multiplication  in  equation  (3 . 8)  and 
separating  the  field  into  its  cartesian  components,  gives  the 
followoing  set  of  three  scalar  equations; 

H^(6,t)  =  1 f  sin^i'sin  a’  +J  cos^p  '  ]  cos  9  dA ' 

X  4-"rC,'’  ji  s  ji 

O  A 

(3.9) 

s  1  ^  ^ 

H  (’S  ,  t)  =  ■: - I  ^[  ( J  -  J„  sin  a  ' )  sin  cos  o'  cos  9 

Y  4TrC,^c»T  i)  s 

o  A 

+  sin  '  sin  9  cosa'jdA'  (3.10) 

H®(9,t)  =  T — [-  sin  9  ( J_sin^'p  '  sin  a'  +J  cos^4)']dA' 
o  A 

(3.11) 

After  substitution  of  dA'  =  p '  d  p'  ds '  ,  and  using  the  odd 
properties  of  some  of  the  terms  appearing  in  the  integral  on 
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: ’  (as  it  has  been  done  previously  in  Appendix  C) ,  the  three 
scalar  equations  representing  the  components  of  the  scattered 
field  reduce  to  the  following  expressions: 


,  t) 

=  cos  V  H^(t) 
o 

(3.12) 

,  t) 

=  -  sin  H^(t) 

o 

(3.13) 

,t) 

=  0 

(3.14) 

3 

where  the  common  term  H^(t)  is  given  by  the  following  integral 

,  S  2-  ,  ~ 

H^(t)  =  p  /  o'ds'  I  4-r(J_  sin  ^j'sin  i’ 

o  0  0  ■  s 

+  J  cos^;')dt'  (3.13) 

Equations  (3.12)  to  (3.15)  represent  the  scattered  field  for 
an  arbitrary  elevation  angle  9.  The  special  case  of  interest 
is  the  baci^scattered  field  where  9  =  t.  For  this  case  there 
is  only  one  component  of  the  scattered  field,  which  is 

(t)  =  -  (t) 

X  o 

and  the  retarded  time  is; 

T  =  t  -  R'/C 


where 


R'  =  z '  . 
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Ail  the  terms  appearing  inside  the  integral  in  equation 
(3.15)  are  independent  of  i'  so  that  the  integration  on  r.  ' 
can  be  carried  out  analytically.  Doing  this  yields  the 
following  expression  for  the  backscattered  field: 


(t) 

s 


-^(J  ,  z  '  ,  r  )  sin  i*  +J  ( 

T  S  ? 


Z 


)  ] t ' ds ' 
(3.16) 


In  equation  (3.16),  the  currents  J  and  J  are  those  given 

S  J 

by  equations  (3.8)  and  (3.9)  respectively,  used  at  time 
-  =  t  -  z'/C  (z*  ^  0)  ,  where  z'/C  is  the  retardation  time  rela¬ 
tive  to  the  origin.  The  backscattered  field  equation  is 
independent  of  i  which  means  that  we  can  integrate  the  time 
derivative  of  the  currents,  along  a  line  on  the  contour  of  the 
body. 


D.  NEAR  ZONE  BACKSCATTERED  FIELD 

In  the  Time  Domain  Scattering  Range  developed  at  NFS, 
the  distance  from  the  antenna  to  the  target  under  test  is  not 
long  enough  to  be  considered  for  all  frequencies  of  interest 
as  a  far  region.  Consequently,  it  is  needed  to  develop  a 
formulation  of  the  scattered  field  in  the  near  zone.  For  this 
case  the  geometry  is  described  in  figure  3.4. 

The  general  equation  for  the  scattered  field  is  given  by: 

^  111  «ri3|1  (3.17) 

A 
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point 


Figure  3.4. 


Geometry  of  the  near- zone 
backscattered  field 


where 


:  =  t  -  R'/C ,  and 

dA '  =  0  '  d  I; '  ds  ' 


r2  = 


+ 


2 


After  carrying  out  the  vectorial  multiplication,  dividing  the 
field  vector  into  its  cartesian  components,  again  dropping 
all  the  odd  terms  in  ip '  and  integrating  analytically  on  t', 
equation  (3.17)  reduces  to: 


(J  R 


-  J  ^  {r  +  z  ' )  1  *3  'ds  ' 
ip  o 


(3.18) 


Equation  (3.18)  is  the  expression  of  the  near  baclcscattered 
field  where: 


R'  =  ( (r  +  z’)^  +  o’^) 

o 

The  currents  J  and  J,(p',z*,t)  are  given  by  equations 

S  l' 

(3.8)  and  (3.9),  being  used  at  the  retard  time  i. 


T  =  t-  (R'-r^)/C  (referred  to  the  origin) 
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IV.  NUMERICAL  SOLUTION  OF  THE  TIME  DOMAIN 
RESPONSE  OF  THE  BODY  OF  REVOLUTION 

In  the  previous  section  the  analytical  equations  for  the 
currents  and  scattered  field  were  derived.  The  equation  ob¬ 
tained  will  have  to  be  solved  numerically,  and  the  first  step 
will  be  to  divide  the  body  into  patches  of  approximately  equal 
curvilinear  surface  area.  This  is  followed  by  the  setup  of 
the  algorithm  for  the  solution. 

A.  DIVIDING  THE  BODY  INTO  PATCHES 

The  geometry  of  the  body  is  completely  specified  by  the 
function  j  (z)  ,  (0  <  z  ^  L) ;  which  represents  the  body  contour 
function  S.  Since  the  integration  is  to  be  carried  over  S, 
the  contour  will  be  divided  into  NS  segments  of  equal  length 
IS.  Consequently,  the  body  contour  function  is  fully  defined 
by  NS+1  points. 

u (  z  ) 

jC 

The  space  samples  of 
segment  and  is  given 

•  — 

1 

z .  = 

1 

The  angle  a  for  each 
at  that  point  and  is 


i  =  1  to  NS+1 

each  segment  cS  is  at  the  middle  of  the 
by 

(,)  +  p  ,)/2 

i  =  i  =  1,  . . . ,NS  (4.1) 

^^1  +  1 

sample  point  is  the  slope  of  the  curve 
approximated  linearly  by 
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a.  =  Arc  tan  (  — 
1  z , 


-)  ,  t  =  i  =  1 ,  .  .  .  ,  NS 


Figure  4.1  shows  the  geometry  of  the  contour  sampling.  Each 
segment  AS  defines  a  ring  at  with  radius  p^.  All  the  rings 
will  be  divided  into  patches  of  approximately  equal  surface 
area,  dA^ ,  according  to  the  following  formulas* 

dA.  =  DS  •  p.  •  A3.  =  AS^ 

1  11 

where  A3^  is  the  angle  between  two  adjacent  sample  points  on 
the  same  ring. 


This  suggests  a  division  of  each  ring  into  an  even  number  of 
patches  NPI  which  will  be 

NPI  =  2xint(^)  (4.3) 


The  function  INT  means  the  integer  value  of  its  argument  so 
that  the  increment  angle  for  the  i*"  ring  becomes 


(4.4) 


Thus,  following  this  procedure,  any  arbitrarily  shaped  body 
of  revolution  can  be  subdivided  into  patches  of  approximately 
equal  surface  area.  Figure  4.2  shows  the  geometry  of  the 
space  sampling  on  the  whole  body. 
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sampling  of  the  whole  body 
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B.  NUMERICAL  SOLUTION  OF  THE  INTEGRAL  EQUATIONS 

The  equations  to  be  solved  are  given  by  Eq.  (3.8)  and 
(3.9)  in  previous  sections. 

For  the  numerical  solution  the  currents  will  be  represented 
by  their  space-time  samples; 

J(,.'>,z,t)  =  J(i,n) 

where  i  is  the  i*"^  space  sample  on  the  contour  function  and 
n  is  the  n^  time  sample 

i  =  1 ,  .  .  .  ,  NS 

n  =  1 , . . . ,NT 

t  =  n.AT 

The  time  being  sampled  at  a  rate  AT.  Similarly  the  incident 
field  is  sampled  and  will  be  denoted  by  its  n-th  time  sample 
at  the  i-th  segment: 

H^(z,t)  =  H^(i,n) 

The  general  procedure  for  the  numerical  procedure  was  already 

mentioned  in  Chapter  II,  and  will  be  reviewed  in  more  detail 

here.  The  incident  field,  H,  is  assumed  to  be  zero  at  all 

times  less  than  some  reference  value  (which  will  be  set  equal 

toi^t).  Hence  the  total  current  on  the  surface  of  the  scatterer 

is  also  zero  prior  to  time  t^. 

At  the  next  time  interval  t,  =  t  +  At  the  incident  field 

1  o 

reaches  the  scatterer  at  z  =  z^^  (first  space  sample)  ,  thus  at 
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time  the  current  will  be  simply  twice  the  tangential  com¬ 
ponent  of  the  incident  field.  As  time  marches  on,  the  incident 
field  covers  more  of  the  body  and  the  current  will  be  due  to 
the  incident  field  plus  the  fields  due  to  currents  at  other 
points  on  the  scatterer  at  earlier  times.  So,  by  marching 
on  in  time,  the  currents  at  each  space  sample,  for  each  time 
increment,  can  be  computed. 

The  currents  appearing  inside  the  integrals  in  equations 
(3.8)  and  (3.9)  are  the  retarded  currents,  and  are  given  by 

J'  (p',2',T)  =  J1  (p ' ,z' ,t-R/C) 

S  f  V  S  ^  V 

The  point  (o',z')  is  represented  by  the  k-th  space  sample 
while  the  time  x  will  be: 

X  =  t  -  R/C  =  n-lt  -  R/C  =  gAT. 

In  the  general  case  x  will  not  be  an  integer  multiple  of 
AT,  so  that  the  current  cannot  be  represented  by  any  of  the 
previously  computed  samples.  In  order  to  find  the  value  of 
the  current  at  time  x ,  there  is  a  need  to  interpolate  between 
the  sampled  known,  discrete  values. 

For  a  better  accuracy  use  was  made  of  the  four  point 
Lagrangian  interpolation  formula  using  two  samples  at  each 
side  of  the  value  gAT  if  t^ ,  t^,  t^,  t^  are  the  four  samples 
used  for  interpolation,  such  that: 

t.  =  t.,  +  At  =  t-  +  2At  =  t,  +  3At 
4  3  2  1 

U  -  ‘34 
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9  3 


t2  =  92 

'^1'  *^2'  '^3'  '^4  integer  values,  and 

=  7  :.t  (g  real  number) 


such  that 


The  interpolation  formula  using  those  four  points  will  be 
given  by ;  [  17 ] 


J 

s , 


(T) 


4  4  T-t. 

n=l  1=1  n  1 
i?^n 


(4.5) 


This  interpolation  formula  must  be  adjusted  in  two  cases: 

1.  Whenever  t^  or  t^  is  greater  than  the  actual  time 
of  computation  t  (in  order  not  to  interpolate  in 
the  future). 

2.  Whenever  t^  or  t^  are  less  than  the  beginning  time 
of  computation  t^  (due  to  causality) . 

The  time  derivative  of  the  currents  will  be  simply  expressed 
as  the  derivative  with  respect  to  x  of  the  interpolated 
value,  which  will  give  the  following  formula: 
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SI(m)  =  sin  [  (m-1)  il3  . 


V  =  Z .  -  2 

i  u 


CJ  =  V  sin  sin  x.  +  k,,  cos  t,  sin  x  .  sina  ,  cosa  . 

K  iK  K  11  K  1 


=  V-CO (m) • sin  a  .  +  ( p  - p . CO (m) ) cos  a. 

1  1  iC  1  1 


F  =  cos  i,  (p  cos(m)  -p.)  +  V  sin  I,  CO  (m) 

Z  K  K  1  K 


NS  represents  the  total  number  of  segments 


NP.  represent  the  total  number  of  patches  at  the 
ith  ring. 


The  total  currents  can  now  be  expressed  in  terms  of  and 


J  ,  which  yields: 


J.(k,n)  =  -2H"(z,  )  +  J  (k,n) 

^  /k  f  li  c? s 


(4.11) 


J.(k,n)  =  -2H  (z,  ,n)  <  sin  (n.)  +J  .{k,n)  (4.12) 

K  K  C  V 

and  are  the  values  of  the  current  components  stored  and 
used  for  computation  of  the  backscattered  far  field,  which 
will  be  expressed  in  numerical  form  as  follows: 


NS 


H"(n)  =  - 


-Q  I  [DJ^  (i,g)  sin  01^ +DJ^  (i  ,g)]lSj^p^  (4.13) 


O  1=1 


where 


g  =  n  - 


z  . 

_ ^ 

clt 


and  DJ^  ^(i,g)  are  given  by  equation  (4.8) 
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Equations  (4.7)  and  (4.8)  used  for  interpolation  of  the 
currents  and  their  time  derivative  are  given  for  the  general 
case  where : 

1.  g^  ■  1  (after  the  incident  field  had  reached  the 
specific  point  i  on  the  scatterer) 

2.  g_^  •-  n  (so  that  no  interpolation  in  the  future  occurs)  . 

In  cases  where  these  two  conditions  are  not  met  the  interpo¬ 
lation  must  be  adjusted  to  include  only  those  points  which 
satisfy  causality  with  no  interpolation  in  the  future. 

The  numerical  formulations  derived  in  this  section  were 
programmed  in  FORTRAN  to  compute  the  currents  and  the  back- 
scattered  field,  due  to  the  excitation  of  a  body  of  revolution 
by  an  impulse  type  field  (or  a  ramp) .  The  input  data  to  the 
program  is  the  function  y(z)  along  with  the  length  of  the 
spatial  segment  IS  and  the  time  increment  tT .  The  program 
is  given  in  Appendix  E. 

C .  NUMERICAL  EXAMPLES 

The  first  example  to  be  considered  is  a  sphere  of  radius 
1  meter.  This  special  example  was  chosen  since  its  result 
appears  in  various  publications  and  comparison  can  hence  be 
made . 

For  the  present  example,  use  is  made  of  the  guidelines  set 
in  Section  0(1. B)  for  the  choice  of  the  incident  field,  the 
sampling  rate  and  the  space  sampling.  The  incident  field  for 
the  case  of  the  impulse  response  is  represented  by  the  gaussian 
shape 
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H 


i 


I  z  ,  t ) 


(  t  - 1 


c 


0 


0 


utherwi se 


where : 


lo" 


which  yields 


IT  =  0.2,  C  sec 


The  pulse  is  described  in  figure  4.3. 


Figure  4.3.  Gaussian  shaped  pulse 
incident  on  a  sphere 

The  pulse  reaches  its  maximum  value  at  the  center  of  the  sphere 


t 


max 


a/C  . 
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The  spatial  segment,  with  15  segments 

IS  =  ^  =  0.209  m  -  C-IT 

la 

The  number  of  time  increments  observed  was  such  that  it 
covers  6  times  the  transit  time  across  the  whole  sphere  (60-AT). 
The  results  obtained  for  the  backscattered  field  are  given  in 
figure  4.4,  on  which  the  results  of  Bennett  [7]  are  superimposed 
for  comparison. 

Observation  of  the  impulse  response  of  the  sphere  shows 
the  following  results: 

1.  The  first  peak  appears  after  a  transit  time  across 
one  radius  of  the  sphere  (which  is  the  time  where  the  inci¬ 
dent  pulse  reaches  its  peak  value) ,  and  is  due  to  the  return 
rrom  the  leading  edge  of  the  sphere. 

The  second  peak  appears  at  about  5.2  a/C  and  is  due  to 
-r.e  r-turn  from  the  back  side  of  the  sphere,  after  a  transit 


--o  rresponds  to  twice  the  free  space  transit  time  up  to 
\r.‘-  snadow  boundary  plus  the  transit  time  around  the  back 
the  sphere.  This  second  peak  is  the  creeping  wave  return. 

..  The  response  as  expected,  starts  decaying  after  some 
time  equivalent  to  5-6  times  the  transit  time  across  the  sphere. 

The  ramp  response  for  the  same  sphere  is  given  in  figure 
4.5.  The  first  positive  incursion  of  the  response  waveform 
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BRCKSCflTTEfiED  FIELD  (RANGEnfl) 


extends  to  about  5  times  the  transit  time  across  the  radius 
which  corresponds  to  1.25  round  trips  across  the  target.  The 
impulse  response  of  a  30  cm  diameter  sphere  was  measured  in 
the  NPS  scattering  range  lab.  The  waveform  obtained  is  shown 
in  figure  4.6.  The  response  appears  here  with  a  reversed 
sign  since  vertical  polarization  is  used  (instead  of  the 
horizontal  polarization  considered  for  the  computations) . 

The  waveform  shows  very  good  agreement  with  the  form  obtained 
for  the  2  meter  diameter  sphere. 

The  second  example  considered  was  the  cone-sphere  whose 
geometry  is  shown  in  figure  4.7.  The  target  contour  was 
divided  into  12  segments  while  all  other  parameters  remained 
the  same  as  for  the  sphere.  The  impulse  response  of  the  cone- 
sphere  target  is  given  in  figure  4.8.  The  waveforms  show  two 
positive  peaks,  corresponding  to  the  leading  edge  and  creeping 
wave  returns.  The  separation  between  them  is  equivalent  to  the 
two  way  free  space  transit  time  up  to  the  shadow  boundary 
plus  the  transit  time  across  the  rear  of  the  sphere. 

[(2’<1.3+-^  <0.7)  =  4.7  light-meters]  The  ramp  response 

of  the  same  cone-sphere  target  is  given  in  figure  4.9,  and 
the  remarks  made  on  the  sphere  apply  here  too. 

The  third  exaunple  was  the  capped  cylinder  target  shown  in 

figure  4.10.  The  contour  of  the  target  was  divided  into  12 

segments,  the  sampling  time  was  AT  =  0.025  m/sec  and 
9  -1 

A  =  4-10  sec  .  The  smoothed  impulse  (gaussian)  response 
of  this  target  is  given  in  figure  4.11.  The  waveform  shows 
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clearly  the  responses  from  the  leading  edge  and  the  rear  of 
the  cylinder,  with  the  same  amplitudes,  while  from  the  flat 
horizontal  region  of  the  cylinder  the  response  is  zero  as 
expected  from  the  analytical  expressions  of  the  field. 

The  ramp  response  for  the  capped  cylinder  is  shown  in 
figure  4.12.  The  first  positive  incursion  of  the  response 
clearly  exhibits  a  shape  comparable  to  the  shape  of  the  target. 
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Radius 


Figure  4.7.  Geometry  of  the  cone-sphere 
target  (half-target) 
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BllCKbCR  i  TERED  FIELD  (RnNCE^hKC) 


0.3 


TIME  (LIGHT-METER) 


Figure  4.9.  Ramp  response  of  the  cone-sphere 


6 
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Figure  4.10.  Geometry  of  the  capped  cylinder 
target  (half-target) 
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BfiCKSCRTTERED  FIELD  (RRNGEhHkC) 
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Figure  4.12.  Ramp  response  of  the  capped 

cylinder  (total  length  '  0.33  m) 


V.  INVERSE  SCATTERING 

A.  PROBLEM  DEFINITION 

Given  the  incident  field  and  the  scattered  field  from  an 
unknown  target,  it  is  desired  in  the  most  general  case  to 
retrieve  the  shape,  size,  orientation,  and,  perhaps,  composi¬ 
tion  of  the  target.  This  is  known  as  inverse  sct.ttering,  or 
target  imaging. 

For  the  present  case,  the  problem  is  restricted  to  a 
metallic  body  of  revolution  illuminated  by  a  field  propagating 
along  the  axis  of  symmetry,  and  the  scattered  field  is  measured 
in  the  backscattered  direction,  as  shown  in  figure  5.1. 

The  solution  of  this  problem  is  approached  here  directly 
in  the  rime-domain. 


B.  FORMULATION  OF  THE  PROBLEM 

In  [6],  Kennaugh  and  Moffatt  derive  the  relationship  between 
the  physical-optics  impulse  response  of  a  body  and  its  area 
function  A(z),  projected  in  the  plane  orthogonal  to  the  direc¬ 
tion  of  propagation  of  the  field  (z-axis) . 

The  impulse  response  is  given  by: 


Fj  (t) 


C  d^A(z) I 
o  dz  ‘  z=ct/2 


(5.1) 


where  r^  is  the  range  from  the  observation  point  to  the 
scatterer,  as  shown  in  figure  5.1.  Integrating  twice  (5.1) 
with  respect  to  time  gives  the  ramp  response  as  follows: 
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fR"=>  "  ‘5.2) 

o  'z=ct/2 

Equations  (5.1)  and  (5.2)  were  proven  to  be  valid  for  bodies 
having  a  monotonic  area  function  with  a  limiting  value  at  the 
shadow  boundary.  This  relationship  between  the  ramp  response 
and  A(2)  is  derived  in  Appendix  F  for  the  body  of  revolution. 

The  total  baclcscattered  field  that  would  be  measured  will 
consist  of  the  physical  optics  term  plus  an  additional  term 
due  to  the  field  induced  interactions  between  the  currents 
at  different  parts  of  the  body.  The  total  current  induced  on 
the  body  is  given  by  equation  (3.1).  The  first  term  of  the 
R.H.S.  of  the  equation  represents  the  physical-optics  current, 
which  we  will  denote  by  Jp ,  and  the  second  term  expressed  by 
the  integral  represents  the  additional  current,  which  we  will 
denote  by  .  The  total  current  will  then  be: 

^total^^'^^  =  Jp(?,t)  t  Jc(r,t)  (5.3) 

Expressed  in  terms  of  (5.3)  the  total  backscattered  field  for 
a  ramp  excitation  will  be: 

-  r4-c  (  Ir'Jp  *  jc>  ^ “r ‘5.4) 

o  A  ^ 

=  H®(t)  +  H®(t) 

p  C 

where : 


H^{t) 


is  the  physical  optics  ramp  response. 
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and 


s 

H^(t)  is  the  additional  response, 

g 

Hj^(t)  is  the  total  ramp  response. 

Using  the  relationship  between  the  physical-optics  response 
and  the  area  function,  {5.4)  can  be  rewritten  as : 

H|(t)  =  j^(z)  +  H^(t)  ,  z  =  ct/2  (5.5) 

which  leads  to 

o^(z)  =  2r^C(H|(t)  -H®(t))  (5.6) 

Equation  (5.6)  is  the  basis  for  the  solution  of  the  inverse 
scattering  problem  for  the  body  of  revolution. 

C.  INVERSE  SCATTERING  SOLUTION  PROCEDURE 

Equation  (5.6)  relates  the  contour  function  of  the  body, 
^  {z) ,  to  its  ramp  response.  Solution  of  that  equation  leads 
to  the  shape  of  the  body.  The  only  term  known  in  equation 
(5.6)  is  the  measured  value  of  the  total  ramp  response  H®(t) 

g 

while  H^(t)  IS  not  known. 

Moreover,  all  the  terms  in  the  equation  are  functions  of 
(z)  ,  hence  there  results  a  non-linear  equation  that  will  be 
solved  by  an  iterative  procedure.  As  a  first  approximation, 
the  unknown  term  H®(t)  will  be  neglected  and  so  the  first 
estimate  of  the  contour  function  will  be  given  by: 

p^(z)  =  [2r^C  H|(t) (5.7) 


71 


Once  a  body  is  defined  its  ramp  response  can  be  computed  since 
we  know  the  incident  field.  The  ramp  response  computed  for 


the  first  estimate  of  the  body  will  be  denoted  by  H_  (t)  and 

^1 

is  expressed  as: 


s 


A  new  estimate  of  the  shape  be  found  using  (5.6) 

and  (5.7) ,  and  is  given  by 


02(2) 


2roC[H“(t) 


H®  (t)  ]  +  phz) 

Rf  1 


(5.9) 


or  equivelently  by 


=  2r^C[2Hp(t)  -  Hj^(t) ] 


:5.10) 


The  new  estimate  serves  to  compute  a  new  ramp  response, 

H  (z) ,  from  which  an  updated  estimate  of  o(z)  is  extracted. 
^2 

Following  this  iterative  procedure,  the  k-th  estimate  of  the 
contour  function  will  be  given  by: 


Pv  (z) 


.k_i(z)  ^  2r^C[H^(t)  -  h| 


:5.ii) 


or  by: 


2  ,  X 

Pj^(z) 


2r  C[k.H®(t)  -  I  h!  (t)l 


;5.12) 


The  iterations  are  continued  up  to  the  point  where  the  differ¬ 
ence  between  the  true  measured  response  and  the  last  computed 
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response  reaches  a  minimum.  Several  runs  of  the  program  for 
various  bodies  has  shown  that,  as  the  iterations  go  on,  the 
mean  squared  error  is  decreasing,  reaches  a  minimum  then  starts 
oscillating.  In  one  case  observed,  the  error  is  diverging 
after  reaching  the  minimum.  Some  representative  results  are 
shown  in  the  next  section. 

The  normalized  mean  squared  error  will  be  defined  as 
follows : 


.2 

'k 


(Hp(t)  -  (t)  ) 

^ 

{Hp(t) 


(5.13: 


2 

£,  is  the  value  that  must  be  minimized  in  order  to  obtain  the 
K 

most  accurate  representation  of  the  shape  of  the  body.  As 
it  is  defined,  the  error  is  equivalent  to  the  difference  be¬ 
tween  the  successive  shapes,  as  shown  by  equation  (5.11).  Once 

the  difference  H„(t)  -  (t)  is  minimized,  the  difference 

^  ‘^k-l 

2  2  ^ 

0j^(2)  -  Cj^  j^(z)  is  also  minimized. 


D.  NUMERICAL  SOLUTION  AND  EXAMPLES 

For  the  numerical  solution  we  will  use  equation  (5.12), 
which  is  more  suitable  for  the  algorithm  developed. 

We  will  incorporate  the  following  terminology: 

Hj^j_j(nAT)  is  the  time  sample  of  the  measured  ramp  response 
while  (nAT)  is  the  time  sample  of  the  computed  ramp  response 
for  the  k-th  estimate  of  the  body,  and  P)^(z)  is  the  k-th 
estimate  of  the  contour  function.  In  discrete  numerical 
form,  equation  (5.12)  becomes: 
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(5.14) 


.  2r^Clk-H^(n4T)  -  V  Hg^(niT)l 


where 


z  =  n  C  IT/l 
n 


and  for  i  =  0, 


2 

"k 


H  (nAT)  =  0.  The  mean  squared  error  is: 
^o 


I  [H  (nAT)  -  H  (nAT)  1^ 
n=l  ^  ^k-1 


n=l 


(5.15) 


where  N  is  the  total  number  of  time  samples  observed  and 
computed. 

The  ramp  response  of  the  body  defir  ad  by  p,  ,  (z  )  will 
be  computed  using  the  integral  equation  solution  to  the  for¬ 
ward  scattering  problem,  as  already  described  in  Chapter  IV. 

The  Fortran  program,  given  in  Appendix  E,  has  been  extended 
with  a  subroutine  named  "SHAPE”,  in  order  to  implement  the 
inverse  scattering  algorithm.  In  this  subroutine  the  shape 
of  the  body  is  computed  iteratively,  using  (5.14). 

The  space  samples  given  by  (5.14)  are  defined  at  equal 
increments  on  the  z-axis,  while  the  program  developed  for  the 
computation  of  the  ramp  response  has  been  proven  to  work  best 
when  the  contour  was  divided  into  equal  length  segments.  Con¬ 
sequently,  one  of  the  tasks  of  subroutine  "shape"  is  to  find 
the  space-samples  at  the  points  delimiting  segments  of 


approximately  equal  length.  This  is  done  using  a  "cut  and 
try"  method  which  converges  to  the  desired  length  of  segment. 

The  iterative  approach  for  solving  the  inversion  equation 
(5.14)  was  tested  in  a  simulative  way  for  the  three  types  of 
bodies  already  considered  for  the  forward  scattering  solution 
in  the  previous  chapter. 

The  ramp  response  of  the  body  is  first  computed  and  is 
used  to  simulate  the  response  that  would  have  been  measured 
(without  noise  or  error) .  Then  the  shape  of  the  body  is  re¬ 
trieved  from  that  ramp  response. 

The  shape  of  the  body  is  retrieved  from  the  first  positive 
incursion  of  the  ramp  response,  which  as  it  will  be  seen 
extends  beyond  the  length  of  the  target  for  the  shapes 
considered . 

The  first  example  was  the  1  meter  radius  sphere.  Figure 
5.2  shows  the  1-st,  3-rd,  6-th  estimates  of  the  contour  func¬ 
tion  of  the  sphere.  The  6-th  estimate  gives  the  least  error 
and  is  shown  superimposed  on  the  original  shape  in  figure  5.3. 

The  second  example  was  the  cone-sphere,  used  in  two  con¬ 
figurations.  The  first  one  when  the  contour  of  the  body  was 
originally  cut  to  12  segments.  Figure  5.4  shows  the  successive 
estimates  (1,3,6)  of  the  contour,  and  figure  5.5  shows  the 
minimum  error  estimate  (estimate  6).  The  second  version, 
when  the  body  was  cut  to  11  segments,  gave  a  much  more  accurate 
result,  as  shown  in  figures  5.6  and  5.7. 
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original  shape 
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Figure  5.3.  6-th  estimate  of  the  contour  of  the  sphere 


oriyinal  shape 
1-st  approximation 
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Figure  5,4.  Contour  estimates  of  a  cone-sphere 
(cut  to  12  segments) 
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Sixth  estima 


The  Last  example  was  the  capped-cylincler .  Figure  5.3 
shows  the  successive  estimates  of  the  shape  (estimates  1,3,4), 
and  figure  5.9  shows  the  4-th  estimate  which  gives  the  minimum 
error . 

For  t.he  case  of  the  capped  cylinder,  after  reaching  the 
minimum  error  at  the  4-th  iteration,  the  imaged  shaped  of 
the  body  started  to  diverge  considerably,  and  appeared  more 
shortened  than  the  actual  shape. 

Common  to  all  examples,  the  error  appears  mostly  at  the 
shadowed  region  of  the  body,  but  in  general  the  shape  is 
retrieved  with  quite  a  good  fidelity,  which  proves  that  the 
iterative  direct  time  domain  solution  is  working,  subject  to 
certain  caveats  which  will  be  discussed  in  the  conclusions. 

E.  INFLUENCE  OF  NOISE 

In  the  examples  considered  in  the  previous  section  the 
only  noise  appearing  in  the  ramp  response  is  the  numerical 
noise  generated  by  sampling  and  computational  round-off.  In 
this  section  the  measurement  noise  will  be  tested  with  a  ramp 
response  to  which  noise  is  added.  In  this  siraulation ,  each 
noise  sample  is  obtained  by  the  sura  of  a  large  number  of  zero- 
mean,  uniformly  distributed  random  num.^^ers,  which,  by  use  of 
tr.e  central  limit  theorem,  approximate  the  statistics  of 
ijissian  noise.  The  noise  samples  are  normalized  in  order  to 
ma.n  within  a  certain  voltage  signal  to  noise  ratio,  and 
added  to  the  computed  ramp  response. 

•  :,o<iy  considered  for  the  noise  test  was  the  sphere, 

.  ;nal  to  noise  ratios  tested  were; 
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Figure  5.8.  Contour  estimates  of  the  capped  cylinder 


in 


V 

20  log  —  =  10,  20,  30  and  40  dB  . 

"noise 

The  sphere  ramp  responses  with  the  noise  added  are  shown  in 
figures  5.10,  5.11,  5.12  and  5.13  for  the  successive  signal 
to  noise  ratios. 

For  the  10  dB  S/N  case,  the  appearance  of  spikes  and  dis¬ 
continuities  on  the  ramp  response  causes  the  shape  of  the 
body  to  diverge  considerably  since  the  first  iteration,  and 
the  shape  cannot  be  recovered. 

For  the  20  dB  S/N  case,  the  computed  shape  starts  to  con¬ 
verge,  with  the  minimum  error  shape  obtained  at  the  third 
iteration.  Then  the  imaged  shape  proceeds  to  diverge  and  to 
shorten . 

Figure  5.14  shows  the  successive  estimates  of  the  shape 
up  to  the  3-rd,  and  figure  5.15  shows  the  sixth  estimate,  which 
is  clearly  different  of  the  actual  shape. 

For  the  30  dB  S/N  case,  the  minimum  error  shape  was  ob¬ 
tained  at  the  5-th  iteration,  as  shown  by  figure  5.16.  Then 
again,  the  computed  shape  diverges  as  shown  in  figure  5.17. 
Common  to  all  cases  of  divergence,  the  wrong  estimates  start 
as  soon  as  the  last  computed  shape  exhibits  a  strong  discon¬ 
tinuity  in  its  contour.  A  very  important  conclusion  from  the 
noise  analysis  is  that  before  starting  the  inverse  scattering 
procedure  the  ramp  response  used  must  be  smoothed  so  that  no 
spikes  or  sudden  discontinuities  appear  on  it. 

We  should  normally  expect  that  the  true  measured  response 
will  have  a  20  dB  s/N  ratio,  so  that  the  problem  of  divergence 
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Figure  5.12.  Ramp  response  of  the  sphere  with 
noise  added  (S/N  =  30  dB) 
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Figure  5.13.  Ramp  response  of  the  sphere  with 
noise  added  (S/N  =  40  dB) 
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must  be  investigated.  In  any  case,  before  diverging  the 
shape  that  gives  the  minimum  error  is  quite  a  good  approxi¬ 
mation  of  the  true  shape. 

An  additional  test  with  40  dB  S/N  was  made,  and  no  signi¬ 
ficant  deviation,  from  the  case  without  noise  added,  was 
observed.  Figure  5.18  shows  the  first  and  6-th  estimates  of 
the  shape.  The  minimum  error  shape  was  the  6-th  estimate. 

Up  to  the  12-th  iteration  the  shape  of  the  body  remained 
stable,  with  minor  variations. 
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VI .  SUMMARY 


A.  SUMMARY  OF  RESULTS 

The  inverse  scattering  method  for  radar  target  imaging 
was  the  main  topic  addressed  in  this  work.  The  inverse 
scattering  algorithm  developed  uses  the  time  domain  solution 
of  the  integral  equation.  The  class  of  targets  considered 
included  the  conducting  bodies  of  revolution  whose  importcint 
symmetry  property  reduces  considerably  the  computation  time 
and  the  computer  storage  needed  for  the  solution  of  the  in¬ 
verse  scattering  problem. 

The  use  of  the  ramp  response  for  retrieving  the  shape  of 
the  target  proved  to  be  a  very  efficient  tool,  and  the  exam¬ 
ples  considered  gave  generally  good  results,  with  initial 
convergence  of  the  imaged  shape.  For  the  case  of  the  long 
cylinder  and  the  sphere  with  low  signal  to  noise  ratios  there 
appeared  an  eventual  divergence  of  the  imaged  shape  after 
many  iterations.  This  is  thought  to  be  due  to  accumulated 
numerical  roundoff  errors  in  repeatedly  solving  the  integral 
equations.  In  general  it  has  been  shown  that  the  use  of 
the  ramp  response  by  itself,  without  any  further  computation, 
gives  a  close  approximation  to  the  shape  of  the  body  but  up 
to  its  shadow  boundary  which  corresponds  to  the  maximum  of  the 
area  function  A{z)  along  the  line  of  sight.  Beyond  that  boundary 
the  body  is,  at  the  first  observation,  elongated  and  the 
iterative  algorithm  is  needed  to  compensate  for  the  errors 
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induced  by  the  wake  or  "creeping  wave"  which  alters  the 
target  response. 

As  the  iterative  computation  process  goes  on,  some  numeri¬ 
cal  instability  that  is  mainly  due  to  a  discontinuity  in  the 
computed  shape  might  arise.  However,  before  diverging,  the 
shape  of  the  body  always  reaches  an  estimate  which  has  a 
minimum  error  relative  to  the  true  shape,  and  is  a  good 
approximation  of  that  shape,  so  that  the  minimum  error  cri¬ 
teria  set  for  stopping  the  iterative  process  appear  to  be  an 
attractive  way  to  overcome  the  instability  phenomena.  As  a 
continuation  of  this  work  it  is  suggested  to  carry  out  a 
detailed  analysis  of  the  error  and  numerical  instabilities 
that  have  been  observed  while  testing  the  iterative  algorithm. 

B .  RECOMMENDATIONS 

The  algorithm  developed  for  target  imaging  by  transient 
inverse  scattering,  has  been  tested  numerically  by  simulating 
additive  physical  noise.  It  must  now  be  evaluated  in  a  real 
laboratory  environment  where  the  true  measured  response  is 
used.  The  trainsient  scattering  range  at  NFS  should  be  used 
tc  test  and  finalize  the  algorithm  for  the  class  of  targets 
addressed  in  this  thesis. 

The  first  step  to  be  done  is  to  obtain  the  measured  ramp 
response  of  the  target.  The  range  uses  an  impulse  generator 
which  exhibits  a  gaussian  pulse  shape  characteristic.  The 
true  ramp  response  can  be  obtained  by  using  either  a  time 
domain  convolution  or  a  frequency  domain  computation  and  I.F.T. 
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The  true  ramp  response  can  be  obtained,  in  principle,  from 
the  gaussian  response  of  the  target  using  the  following 
relationship ; 


where ; 


Hr(.) 


G  {‘i) ) 


2 


HrIoi)  is  the  Fourier  transform  of  the  true  ramp 
response  of  the  target, 

H  (oj)  is  the  gaussian  impulse  response  of  the 
target, 

G(u))  is  the  Fourier  trcinsform  of  the  transmitted 
gaussian  impulse, 

and  the  ramp  response  will  be; 

HR(t)  =  f‘^[Hr(w)]  . 

Very  special  care  must  be  given  to  the  choice  of  space 
sampling  and  time  sampling  intervals.  The  time  sampling  is 
dictated  by  the  sampling  rate  of  the  oscilloscope,  and  the 
space  Scunpling  must  be  done  at  a  high  enough  rate  such  that  in 
every  case  the  distance  R  between  the  closest  spatial  samples 
is  related  to  the  time  seimpling  T  by: 

R/C  ^  AT  . 

The  next  step  is  to  translate  the  FORTRAN  program  given  in 
Appendix  E  to  fit  the  computer  of  the  lab.  There  is  a  need 
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to  reduce  to  a  minimum  the  computation  time  since  the 
laboratory  microcomputer  is  very  slow  compared  to  a  large 
mainframe.  The  use  of  some  additional  symmetries  of  the 
body  can  certainly  help  in  this  task. 

Finally,  it  is  necessary  to  minimize  the  noise  in  measure¬ 
ments,  and  smooth  the  response  to  be  used  for  shape  computation. 

C.  AREAS  FOR  FURTHER  STUDY 

The  ultimate  objective  of  this  research  progrcim  is  to 
aid  in  determining  the  feasibility  of  developing  future  radar 
systems  having  the  innate  ability  to  classify  targets  using 
their  transient  scattering  response. 

In  the  most  general  case  the  target  will  be  of  arbitrary 
shape  cind  the  use  of  the  iterative  method  described  in  this 
thesis  would  be  highly  impractical  and  prohibitive  for  com¬ 
plex  targets,  due  to  two  principal  factors.  The  first  is  the 
time  needed  to  compute  the  response  of  a  complex  body,  and 
the  second  is  that  the  area  function  obtained  from  the  ramp 
response  is  not  sufficient  to  determine  in  what  way  the 
response  is  to  be  computed. 

The  best  technique  to  deal  with  such  targets  would  involve 
a  way  that  does  not  require  the  computation  of  the  integral 
equation,  and  instead  will  use  only  the  measured  impulse 
response  {or  ramp  response)  to  discriminate  between  targets. 

The  first  solution  involves  measuring  the  ramp  response  at 
multiple  look  auigles.  The  response  waveform  will  be  used  up 
to  the  shadow  boundary  to  find  the  area  function  along  each 
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line  of  sight.  The  target  image  can  be  approximately  recon¬ 
structed  from  the  projections  of  the  area  functions.  Das 
and  Boerner  [18],  suggested  a  very  interesting  approach  for 
target  shape  estimation  using  these  projections.  This  method 
seems  attractive  and  is  felt  that  it  could  lead  to  a  good 
solution  for  target  imaging.  Its  main  disadvantage  is  the 
need  for  multiple  look  angles,  requiring  either  interfaced 
multiple  radars,  or  long-time  dwelling  with  one  radar. 

An  alternative  approach  is  to  classify  targets  via  the 
complex  poles  and  residues  of  their  time-domain  signatures. 

Once  the  poles  and  residues  are  extracted  they  could  be  com¬ 
pared  with  stored  library  data  corresponding  to  known  targets, 
and  thus  would  provide  a  classification  of  the  target  under 
test.  This  catalogue  approach  is  also  being  pursued  using 
the  NFS  transient  range,  as  well  as  through  other  independent 
efforts,  and  may  prove  to  be  the  most  practically  implementable 
technique  available. 

Another  area  which  could  form  the  subject  of  future  re¬ 
search  includes  the  time  correlation  of  the  shape  of  the 
target  with  its  motion. 
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APPENDIX  A 


DERIVATION  OF  THE  SCALAR  INTEGRAL  EQUATIONS 

In  Chapter  III,  equations  (3. 3)  and  (3. 4)  involve  some  vector 
multiplications.  The  first  appearing  in  the  integrals  is: 


a  X  J '  X  R  = 
n 

(a  X  a  '  : 
n  s 

X  R)  .  j  '  +  (a^  X  a  ' 
s  n  P 

X  R)  j; 

•p 

(Al) 

We  note  the  following  identities: 

a  X  a  '  X  R 
n  s 

=  a' (a  • 
s  n 

R)  -  R(a  -a') 
n  s 

{A2) 

a  X  a  '  X  R 

n 

“  al(a_* 

n 

R)  -  R(a  .aM 
n  p 

(A3) 

Expressing  R,  a 

, '  and  a ! 
s  P 

in  terms  of  the 

coordinates 

system  gives ; 

R  = 

+  a  (R*a  ) 
s  s 

- 

a  (a* -a  ) 

P  P  P 

/s 

+  a  (a ' -a  ) 

S  (p  s 

(A4) 

"s  = 

a  .  (a  ' -a  ) 

S  S 

y\ 

+  a  (a' -a  ) 
s  s  s 

For  the  body  of 

revolutions  we  note  the 

following 

relationships : 

‘P  (p 

=  cos  0 

^  /\ 

=  -  sin 

a -sin  3 

a'  -a 

Ip  n 

=  -  cos 

a* sin  3 

a '  •  a 
s  s 

=  sin  a 

sin  a 'cos  8  +  cos  a  cos  a 

1 
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sin  a*  sin  6 


X  a  . 

s  t> 


a '  -a 
s  n 


=  cos  a  sin  a 'cos  6  -  sin  a  cos  a' 


where 


6  =  4)’  - 


The  vectors  R,  a  ,  a  ,  a,  are  given  in  cartesian  coordinates 
n '  s  t) 


R  =  a  (p  cos  o-L)'  cos>')  +  a(psin4)-p'  sin  4)') 

X  y 


-  (z '  -  z) 


a  =  a  (cos  a  -cos  4))  +  a  (cos  a  -sin  (p)  -  a  sin  a 

n  X  y  z 

a  =  a  (sin  a  -cos  4))  +  a  (sin  a  ‘Sin  4> )  +  a_  cos  a 

s  X  y  z 


a^  =  -a.^(sin  a)  +  a^^cos  (p 


Carrying  out  the  dot  products  gives: 


R-a  =  cosa(p-p'cos3)  +  (z'-z)sina 
n 


R-a  =  sin  a(p  -o'cos  3)  -  (z'  -z)cos  a 
s 


p '  sin  6 


Using  the  relationships  given  by  the  sets  of  equations  (A4) 


and  (AS),  and  simplifying,  equations  (A2)  and  (A3)  become: 
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a  a  '  X  R 
n  s 


=  a^lsin  [  ( z  ' -2)  sinasina  '  +  ocosotsina  ' 


'  cosa  '  sina]  }  +  a^  •'  ( p-p  'cosB)  cosa  ' 


+  { z ' -z) cosBsina ' } 


a  xa'xR  =  a^  [  (z '-z)  sinacos  B+ cosa  (pcos8-p  ' )  ] 

lip  (P 


+  a  [-sinB ( z ' -z) ] 
s 


Defining 


V  =  z '  -  z 


=  (z  '  -  z)  cos  6s  in  a'  +  cos  a'(p-p  'cos  6 ) 


U  =  (z  '  -  z)  sin  a  sin  a  '  +  p  cos  a  sin  ct '  -  p  '  cos  a  '  sin  a 


=  cos  a (pcos  6 -p ' )  +  sin  a  cos  6 (z '  -  z) 


produces  from  (A6)  and  (A7) 


a„  X  a '  X R  =  a,  u  sin  6  +  a^-F, 
ns  si 


a„  X  a '  X R  =  a.  F_  -  V-sin  3 
n  p  <p  2  s 


Note  that  V,  F^,  F^  are  even  functions  of  8  and  U  is  inde¬ 
pendent  of  S.  Next  we  have  to  carry  the  second  vector  multi¬ 
plication  ; 


"  — i  ^  i  .  ^  i 

a  X  H  =  -a  H  sin  <{)  -  a .  H  cos  i)  sin  a 

ns  <p 


(AlO) 
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where  H  is  given  in  cartesian  coordinates  byi 


After  substituting  equations  (A8) ,  (A9),  and  (AlO)  into 

equations  (3.3)  and  (3.4)  we  obtain  the  scalar  expressions 
for  the  components  of  the  current  in  the  s  and  -p  directions 
as  follows: 


Jg  =  -2H^sin  $  +  i  |:pl  [J^F^  -  Vsin  6]dA'  (All) 

A  R 

=  -2H^cos  psin  a  +  i  1^] 

A  Xx 

• [Jl  U  sin  6  +  J'F_]dA' 
s  'P  ^ 


(A12) 


APPENDIX  B 


DEMONSTRATION  OF  THE  SINUSOIDAL  VARIATION 


OF  THE  CURRENTS  WITH 


Maxwell's  first  equation: 


7  <  H  =  E  +  J 


(Bl) 


In  the  case  of  interest,  due  to  symmetry,  the  field  H  is 
independent  of  J  and  is  given  by 


H  =  a^  H^(p,z,t) 


or  in  cylindrical  coordinates 


H  =  H  cos  0  a  -  H  sin  ■!> 

o  o  '•? 


3H. 


■p  6Z 

9H 

o 

-  a  sin  i  '-T — 
z  3p 


(B2) 


The  currents  on  the  body  expressed  in  cylindrical  coordinates 
are : 


J  =  a  sinaJ  +  a,J,+a  cos  a  J 
P  S  >  ({)  z  s 


The  electric  field: 


E  =  a  E  =  a  sin  4)E  +  a^  cos  <p  E 

yo  p  od)  '^o 


(B3) 


(B4) 


Substituting  equations  (B2)  ,  (B3)  and  (B4)  in  Maxwell's 
equation  (Bl)  aind  equating  the  terms  in  the  same  direction  gives 
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sin  i-J  =  [—  4— )  -  E  ]sin  p 

s  dz  o  o 

The  term  inside  the  brackets  is  independent  of  p  so  that  we 
can  write 


sin  o(  J  = 


(t;,z,t)  sin  -p 

°1 


Similarly 


3H 

o  „  , 

—  -  E  )  cos  p 
3  2  o 


J  =  J  (p,z,t)  cos  p 

P 


and  finally 


3H 

J  cos  a  =  -  sin  a  ::: - 

S  .Cq 


(B5) 


(B6) 


=  J  (p,z,t)  sin  p  (B7) 

Equations  (B5) ,  (B6) ,  (B7)  show  clearly  that  the  currents 

Jg  and  vary  sinusoidally  with  the  angle  p  so  that  we  can 
generalize : 


Jg(p,z,t,p)  =  Jg(p,z,t)  sin  p 


J  p  ( p  ,  z  ,  t ,  P  ) 


J  ^  ( p , z , t)  cos  p 
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APPENDIX  C 


SIMPLIFICATION  OF  THE  SCALAR  INTEGRAL  EQUATIONS 

In  Chapter  III  the  scalar  equations  for  the  currents 
components  were  given  by  equations  (3,4)  and  (3.5),  In  these 
equations  we  had: 

dA  '  =  c '  ds '  dB 


and 


Js  =  J  sin  •:> 
s 

J'J  =  J  cos  5 
? 


Substituting 


J  sin  P  = 
s 


in  equation 
-2H^sin  p  + 


(3.4) 


gives 

s 

'ds  ' 


0 


2-r  ,  ,  . 

1  r  1 

^2^C  3t 


•  [ J  '  sin  p'F,  -  J  '  sin  S  V  cos  P  '  ]d3  (Cl) 

s  IP 

As  noted  in  Appendix  A,  the  expressions  V,  F^  and  F^  appearing 
in  the  integral  are  even  and  periodic  functions  of  c.  We 
recall  that 

R  =  (p^  +  c  '  ^  -  2POCOS  3  +  (z  '  -  z)  ^) 


and 


T  =  t  -  R/C 
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so  that  R,  Jg(p',z,T)  and  J'(p*,z’,t)  are  also  even  and 
periodic  functions  of  Finally  U  is  independent  of  3. 

We  will  concentrate  on  the  integral  over  t  in  order  to 
reduce  it  to  its  simplest  expression.  Since 

4)'  =  '!)  +  S 


we  can  write 


sin  d'  =  sin  j  cos  B  +  sin  B  cos  0 
cos  '  =  cos  cos  3  -  sin  •))  sin  6 


and  substitute  this  into  the  integral  over  3  in 
equation  (Cl).  The  integral  over  S  for  a  complete  period,  of 
all  terms  that  contain  sin  S  will  be  zero  since  those  terms 
are  periodic  and  odd  functions  of  3.  All  the  other  terms  will 
produce  a  non-zero  integral  and  will  be  retained. 

After  this  simplification,  equation  (Cl)  will  become: 


J  sin  5 
s 


i.  1  ^ 

-2H  sin  <ti  +  sin  (p  /  o'ds' 

2tt  ,  ,  .  .  - 

•  I  ^  ^  +  [J'F-  cos  8  +  VJ’sin^SldS  (C2) 

n  OT  R  si  P 


which  after  repressing  the  common  term  sin  p,  leads  to  the 
final  equation; 


2tt  . 


J  =  -2H^  +  ^  J  p'ds'  / 

s  2Tr  r2  C  ..T  R 


•  cos  8  +  VJ^sin^6]d8  . 


(C3) 
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Following  the  same  reasoning  for  equation  (3.5),  and  repressing 
the  common  term  cos  i  leads  to  the  following  final  equation: 


1  ,""l,l 


=  -2H^  Sin  a  +  /p’ds*  / 


0  R 


2 

[J'  U  sin  3  +J'F_cosbld3 
s  i  2 


In  both  equations  (C4)  and  (C3)  the  currents  are: 


S  ,  y  a  /  P 


T  being  the  retardated  time  t  -  R/C,  and  the  incident  field 


=  H^(2,t) 
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SYSIN  DO  DATA 


APPENDIX  D 


4c 

*  TRANSIENT  RESPONSE  OF  A  STRAIGHT  THIN  WIRE 

4c 

«  4t «  «  «  «  «  «  «»  4c  *  «  4c «  «  4c  *  «  4c  4c «  «  *  «  »  »  »  «  «  4c «  *  *  4c «  4c «  4c «  «  4c  4c «  *  »  «  *  4c  4c 


P'JRPQSE: 

4i4>4c4<*4c4i  4c 

THIS  PROGRAM  COMPUTES  THE  TRANSIENT  TIME  DOMAIN 
RESPONSE  OF  A  STRAIGHT  THIN  WIRE. 

THE  PROGRAM  COMPUTES  THE  CURRENTS  INDUCED  ON  THE 
WIRE  dY  AN  INCIDENT  IMPULSE  FIELD, AND  THE  SCATTERED 
FIELD,  FOR  ANY  ANGLES  OF  INCIDENCE  AND  SCATTERING. 


USAGE: 

4c«*4c*4c 

THE  PROGRAM 
OF  THE  INITIAL 
OF  THE  DESIRED 


CAN  BE  USED  AFTER 
PARAMETERS  OF  THE 


PROPER  DEFINITION 
MIRE, AND  DEFINITION 


ANGLES  OF  INCIDENCE  AND  REFRACTION. 


THE  WIRE  IS  ALIGNED  WITH  THE  Z-AXIS  STARTING  AT  Z«0 
THE  ANGLES  ARE  MEASURED  WITH  RESPECT  TO  THE  2-AXIS 


PARAMETER  DEFINITION 

4c  4c4c4c  *4c  4c4c  *4c4c*  4c  4c«4c4c4c 

dETA{A00,40)  : ARRAY  DEFINING  THE  CURRENTS  FOR  UP 

TD  40  SEGMENTSAND  400  TIME  INCREMENTS. 

GAMMA! 400, 401  : ARRAY  DEFINING  THECHARGES  FOR  UP 

TD  40  SEGMENTS  AND  400  TIME  INCREMENTS. 

NS:  NUMBER  OF  SEGMENTS  IN  WHICH  THE  WIRE  IS  DIVIDED 

SL:  TOTAL  LENGTH  OF  THE  WIRE  IN  METERS. 

A:  RADIUS  OF  THE  WIRE  IN  METERS. 

ITS:  NUMBER  DF  TIME  INCREMENTS. 

ITS  MUST  BE  A  MULTIPLE  OF  80  FQR  PLOTTING 
CONVENIENCE 

TI:  ANGLE  OF  INCIDENCE  IN  DEGREES. 

TS:  ANGLE  OF  REFRACTION  IN  DEGREES. 

DZ:  LENGTH  OF  A  SEGMENT  (L/NS) 

DT:  TIME  INCREMENT  IN  SECONDS  (OT*OZ/C) 

AA:  GAUSSIAN  IMPULSE  PARAMETER. IN  SEC-l. 
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t:-iax;  time  delay  of  the  impulse  from  the  first  time 

CF  INCIDENCE  TO  ITS  PEAK. ( TYPICALLY  DTONS/ZI 

El;  INCIDENT  FIELD  (GAUSSIAN  PULSEK 

ES(40D1;  ARRAY  CONTAINING  THE  COMPUTED  SCATTERED 
FIELD  UP  TO  400  TIME  INCREMENTS. 

DELAZ  :  TIME  DERIVATIVE  DF  THE  MAGNETIC  POTENTIAL. 

DELPHI;  TIME  DERIVATIVE  DF  THE  ELECTRIC  POTENTIAL. 

C;  SPEED  OF  LIGHT. 

Cl:  CONSTANT  (MU/4*PII. 

C2;  CONSTANT  (  1 /4*PI  ♦EPS  )  . 

X(30I,Y(30):C0aRCINATES  REQUIRED  FOR  SUBROUTINE 
• PLOTP* . 


REMARKS; 

THE  ANGLES  OF  INCIDENCE  AND  REFRACTION  DESIRED  ARE 
ENTERED  IN  DATA  CARDS  ACCORDING  TO  FORMAT  2F6.1. 

EACH  DATA  CARD  CONTAIN  ONE  ANGLE  OF  INCIDENCE 
AND  ONE  ANGLE  OF  REFRACTION. 

THE  COMPUTATION  IS  DONE  FOR  AS  MAIJY  DATA  IMPUTS  AS 
there  ARE  DATA  CAROS. THE  LAST  DATA  CARD  MUST  CONTAIN 
THE  DATA  NUMBER  l.O  IN  ORDER  TO  TERMINATE  THE  PROGRAM. 


SUBROUTINES  USED: 

««♦♦♦♦♦♦♦♦♦♦♦♦♦♦ 

THE  PROGRAM  USES  A  LIBRARY  SUBROUTINE  NAMED  PLOTP 
FOR  THE  PLOT  OF  CURRENTS  AND  FIELD  AS  A  FUNCTION 
OF  TIME. 


METHOD  OF  SOLUTION: 

ft AAA  Aft  AAA A  AA AAA AAA 

THIS  IS  the  SOLUTION  OF  THE  E.F.I.E.  EQUATION 
APPLIED  TO  A  THIN  W IRE , ACCORD  I NG  TO  THE  METHOD 
DESCRIBED  BY  HARRINGTON  AND  SAYRE  IN  *TIME  DOMAIN 
RADIATION  ANO  SCATTERING  BY  THIN  WIRES  •  I  IN  'APPL. 
SCI. RES.  26*  SEPT  19721 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 


c 

c 

c 

c 

c 

c 


c 

c 


c 

c 

c 

c 


c 

10 

8 

c 

c 

9510 

c 

c 

c 

9601 


c 

c 

c 


c 

c 


♦  * 
*  MAIN  PROGRAM  * 

««****«*  «4i«4i4<4iti*4< 


DIMENSION  BETA!  400 t40 It GAMMA ( 400 «40  If  ESI  400) 
DIMENSION  X(400) ,Y(400I 


CONSTANTS  INITIALIZATION 

C=3.0E+8 

CI=l.0E-7 

C2=<).0E+9 

PI  =  3.  141592654 


PARAMETERS  INITIALIZATION 

COMPUTATION  PARAMETERS 

NS=40 

ITS=400 

NSl*NS-l 

WIRE  PARAMETERS 


A=7.75  E-4 

DZsSL/NS 

DT«02/C 


GAUSSIAN  IMPULSE  PARAMETERS 

AA*0.57  E*9 
TMAX«0T*18.0 

CONTINUE 

READING  OF  DESIRED  ANGLE  OF  INCIDENCE  TI 
AND  ANGLE  OF  REFRACTION  TS. 

REA0(5, 95101  TI,TS 
FORMAT! 2F 6. II 
IFITI.EO.l.OI  GO  TO  999 

WRITE  HEADING 

WRITE(6, 96011  TI,TS 

FORMAT! IHlflOX, 'CURRENTS  INDUCED  WHEN  FIELD  INCIDENT', 
2'AT  ANGLE' f IF6.1,// flOXt'FAR  FIELD  SCATTERED  TOWARDS', 
3'ANGLE' ,IF6. I,////, 8X,' TIME', 15X,' CURRENT  ! AMP  I ',///) 

TRANSLATION  FROM  DEGREES  TO  RADIANS, 

TIR»PI*TI/180.0 

TSR«PI*TS/180.0 

AI»GOS(TIRI 

AR«COS!TSR) 


112 


C  START  OF  CURRENTS  COMPUTATIOr4 

C 

c 

Sl*0.5 
S2=-0.5 
S3=< A/DZ»**2 

FO»ALaG((  S1*SQRT(  Sl**2^S31l/IS2*SQRTIS2**2  +  S3))  ) 

C 

C  INITIAL  CONDITIONS  SETTING 

C 

DO  100  L=1,NSL 
Ll*L-l 

IF(TI.LE.90.0I  30  TO  120 
L2=L-NS 
3J  TO  130 
120  CONTINUE 
L2=L 

130  CONTINUE 

El*  EXP((-AA**21*(0T*( 1-L2*AI »-TMAX)**2)  ♦ S I N (T IR) *3 . 0 
B£TA(1,L»  »£I/(2*C1*FQ/DT» 

IFIL.EQ.ll  GO  TO  110 

GAMMA! l,L i*-( 1 .0/C ) *1 3  ETA ( I ,L I -SET A( 1 ♦ L 1 ) I 
GO  TO  100 
C 

C  BOUNDARY  CONDITIONS  SETTING 

C 

110  CONTINUE 

GAMMAd.l  >=-BETA  (1,1  I/C 
100  CONTINUE 

GAMMAdfNSI*  3ETA  (1,NS1)/C 
C 

C  START  ALL  CURRENTS  COMPUTATION 

C  AT  EACH  TIME  INTERVAL, AND  FOR  ALL  SEGMENTS 

C 

DO  200  J*2,ITS 
C 

DO  300  L=1,NS1 
C 

DELPHI  =0.0 

0ELA2  *0,0 

DO  310  K*1,NS1 
K1=IA8S{L-K) 
i^2*J-Kl 
Sl*Kl«-0.5 
S2*Kl-0.5 
S3=(A/0ZI **Z 

F  »ALOGd  S1*-SQRT(  Sl**2«-  S3)  I  /(  S2«-SQRT(  S2**2«-S3  ) )  I 
IF(Kl.EQ.O)  GO  TO  320 
IF(K2.LE.0J  GO  TO  310 
OSLAZ  *OELAZ  ♦8ETA( K2, K) ♦F 

320  CONTINUE 
K3»K2-l 

IF(K3.LE.0I  GO  TO  310 
OELAZ  »OELAZ  -BETA(K3, K ) *F 

310  CONTINUE 

OELAZ  »OELAZ  *Cl 

DO  330  K»l,NS 
X2*IASS  d.-KI 
K1*IABS  (L-K4-1) 

K3» J-l-Kl 
K4» J-1-K2 

IF(K3.LE.0I  GO  TO  340 

Sl*Kl+0.5 

S2»Kl-0.5 

F  *ALOGI  (  SI  ♦•SORT!  SI  ♦♦2+ S3)  I /I  S2^SQRTI  S2**2*S3 ) )  ) 

DELPHI  -DELPHI  ♦GAMMA ( K3 , K) *F 

C 
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340 


330 

C 

C 

c 

c 

c 

c 

331 

332 

333 

334 


350 

300 


C 

C 

c 

c 

c 

9690 

200 

C 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 


c 

c 


CONTINUE 

IF(K4.LE.0)  go  to  330 

Sl=K2^0.5 

S2=K2-0.5 

F  aALJGK  S1«-SQRT(  Sl**2<-S3)l/(S2+SQRT{S2**2  +  S3)n 
■OELPHI  ^DELPHI  -GAMMA IK4t  K )  ♦F 

CONTINUE 

OELPHI  aOELPHI  *€2 

IF(T1,LE.90.0I  30  TO  331 

FOR  ANGLES  OF  INCIDENCE  MORE  THAN  90  DEGREES 


L2=L-NS 
GO  TO  332 

FOR  ANGLES  OF  INCIDENCE  LESS  THAN  90  DEGREES 

CONTINUE 
L2  =  L 

CONTINUE 

IF  (J.GT.50  )  GO  TO  333 

EI=  £XP(  (-AA#>i‘2l*(  J*0T*<l-L2*An-TMAX)**2»  *S1N(TIR)*3 
GO  TO  334 
CONTINUE 
EI*0.0 
CONTINUE 

BETA(J,L)=( EI-DELAZ 
Jl=J-l 
L1=L-1 

IFIL.GT.I I  GO  TO  350 
GAMMA! J»L)«  GAMMA! JltL 
GO  TO  300 
CONTINUE 

GAMMA! JfL)=  GAMMA! JltL 
CONTINUE 

GAMMA! J I NS J*  GAMM A ! Jl ,NS ) ♦BET A! J t NS H /C 
KK*NS/2 


/DT-OELPHI 

»-B£TA!J,L 

)-BETA!J,L 


/OZ)/!Cl*FO/OT) 

)/C 

)/C  ♦BETA! J,Ll)/C 


rtRITE  THE  CURRENT  INDUCED  AT  THE  CENTER  OF  THE  WIRE 


WRI TE(6t9690l  JfBETA!J,KK» 
FORMATIlOX, I4f21X, E12.5) 

CONTINUE 

END  OF  CURRENTS  COMPUTATION 


START  OF  SCATTERED  FIELD  COMPUTATION 


DO  400  LP»lfITS 
ES1LPI*0. 0 
DO  500  K*lfNSl 
IF  ITS. GT. 90.01  GO  TO  501 
IF  ITS. EQ. 90. 01  30  TO  503 

FOR  REFRACTION  ANGLE  LESS  THAN  90  DEGREES 

KP«K 

N3-LP-1-! KP-O.SIOAR 
NQP*LP-1- !KP^0,5)*AR 
GO  TO  502 
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c 

c 

501 


c 

c 

c 

503 


502 


510 


520 


530 

500 


400 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


600 

9630 

9600 

c 


FOR  KtFRACTION  ANGLE  MORE  THAN  90  DEGREES 

CONTINUE 

KP=K-NS 

NQP=LP-1-(KP-0.5)*AR 
NQ=LP-L-(  KP-t-O.SlOAR 
GO  TO  502 

FOR  REFRACTION  ANGLE  EQ.  90  DEGREES. 

CONTINUE 

NQ=LP-l 

NQP*NQ 

CONTINUE 

NOl*Ng-l 

NQ2=  NQP-1 

IF(NO.LE.O)  GO  TO  500 
IFINQ.EQ. 11  GO  TO  510 
IF( NQ.EU.NQP)  GO  TO  520 
Z3=(LP-l-N0)/AR  *-0.5 
EPS=0Z*<ZQ-KP1 
QEL^Dl-EPS 

ES(LP)*ES(LP) ♦EPS*(BETA(NQ,K)-BETA(NQ1,K) ) 
IF(NQ.EQ.2)  GO  TO  530 

ES<LP»  =  ES(LP» ♦OEL*  ( BET  AT NQP,  K» -BETA(NQ2 ,  K >  » 

GO  TO  500 

CONTINUE 

ES(LP)*ES(LP»  BETATNQtK) 

GO  70  500 
CONTINUE 

ES<LP)*ES  (LPI+OZ^TBETATNQfKI-BETAINQUK)  ) 

GO  TO  500 
CONTINUE 

ES<LPJ»ES(LP)+0Z*2*  BETATNQPtK) 

CONTINUE 

CT=-SIN(TSR>*(1.0  E-7I/0T 
eSTLP)=CT*ES( LP) 

CONTINUE 


END  OF  FIELD  COMPUTATION 


PLOT  OF  CURRENT  AT  CENTER  OF  THE  WIRE 
PLOT  EVERY  MM  POINT. 

MM=ITS/100 
00  bOO  I=»MM,ITSfMM 
K=I/MM 
X(KJ  «I*OZ 
KK=NS/2 

VTK)*  BETATI.KK)  *1000.0 
CONTINUE 
.MM=  ITS/8  0 
WRITEI6,9630) 

FORMATT  IHl) 

.M00CUR»0 

CALL  PLOTP<X,Y,80  ,MOOCUR| 

WRITE  (6,96001  TI 

FORMATT// //,5X, 'CURRENT  INDUCED  AT  CENTER  OF  WIRE',' 

2  SCA^TERER*,  /,5X, 'FIELD  INCIDENT  AT  ANGLE* , 1F6. 1,///, 

3  lOX, 'X-AXIS  TIME  IN  LIGHT-METERS', 

4  //,10X,* Y-AXIS  CURRENT  IN  MILLIAMPS') 
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PLOT  OF  SCATTERED  FIELD 


9640 


9620 


Oa  700  I  =  MM, ITS, MM 
K=I /MM 

Y(K)»ES(I)  ♦1000.0 

CONTINUE 

WRITE  (6,9640) 

FORMAT(Hl) 

MOOCUR=0 

Call  plotp(x,y,30  ,mqocur) 

WRITE  (6,9620)  TS 

FORMAT!////, 5X, 'NORMALIZED  FAR  FIELD  BACK  SCATTERED', * 
2BY  WIRE',  /,5X,*AT  ANGLE* , 1F6 . 1 , / / , 1  OX, • X- AX  I S  TIME '  , ' 
3  LIGHT  METERS' .// ,10X, 'Y-AXIS  FAR  FIELD  MILLIVOLTS') 


END  OF  PROGRAM  GO  BACK  TO  NEXT  ANGLE  IF  ANY 


GO  TO  10 
CONTINUE 
STOP 
Ei^^0 


EXAMPLE  OF  DATA  CAROS 


90.0 


nooorioooor>noooooonr>onr)nnnoor)Or>r>ooooooooo 


C 


APPENDIX  6 


*  « 

♦  TIME  DOMAIN  RESPONSE  AND  INVERSE  SCATTERING  ♦ 

*  <■ 

#  FOR  A  PERFECT  CONDUCTOR  dOOY  OF  REVOLUTION  * 

*  * 


***********************  ***********i^************* 


NAME  OF  PROGRAM  :  GAIL 

**************** 


PURPOSE: 

******** 


THIS  PROGRAM  COMPUTES  THE  TRANSIENT  TIME  DOMAIN 
RESPONSE  OF  A  PERFECT  CONDUCTOR  BODY  OF  REVOLUTIONt 
AND  applies  IT  TO  THE  SOLUTION  OF  THE  INVERSE 
SCATTERING  PROBLEM. (COMPUTATION  OF  THE  SHAPE  OF  THE 
BODY  FROM  ITS  TIME  DOMAIN  RAMP  RESPONSE) 


117 


onoonoooooooooonoooooooooonooriorioonoooonnnortnorvoonnooooooodoooooonn 


USAGE: 


THE  PROGRAM  CAN  8c  USED  IN  ONE  OF  FOUR  WAYS 
DEPENDING  OF  THE  DESIRED  OUTPUT  AND  CONSEQUENTLY  OF 
THE  VALUE  ASSIGNED  TO  THE  INITIATING  PARAMETERS  'INC* 
AND  ’INVER*. 

THE  POSSIBLE  USES  ARE: 

1.  COMPUTATION  OF  THE  TRANSIENT  TIME  DOMAIN  IMPULSE 
RESPONSE  (BACKSCATTEREO  FIELD  ),WHEN  THE  INCIDENT 
FIELD  IS  A  GAUSSIAN  IMPULSE  .(INC=1) 

THE  CURRENTS  INDUCED  ON  THE  BODY  CAN  BE  MADE  AVAILABLE 
AT  THE  OUTPUT  IF  DESIRED  (FOR  EMP  PROBLEMS). 

THE  INPUTS  REQUIRED  ARE: 

A.  THE  SPACE  SAMPLES  OF  THE  CONTOUR  FUNCTION  OF 
THE  BOOY,fia(Z),ENTERO  IN  OATA  CARDS. 

B.  THE  NUMBER  OF  SEGMENTS  SUBDIVIDING  THE  CONTOUR 
OF  THE  BODY  (NS<20). 

C.  THE  LENGTH  OF  A  TIME  INCREMENT  IDT). 

D.  THE  NUMBER  OF  TIME  SAMPLES  NEEOED  I  ITS<*  60). 

E.  THE  length  OF  EACH  SEGMENT  0S(20). 

F.  THE  PARAMETERS  DEFINING  THE  GAUSSIAN  IMPULSE. 

2.  COMPUTATION  OF  THE  RAMP  RESPONSE  OF  THE  BODY 
WHEN  THE  INCIDENT  FIELD  IS  A  RAMP.(INC»2  AND  INVER* I ) 

THE  INPUTS  REQUIRED  ARE: 

SAME  AS  A  THRU  E  LISTED  ABOVE. 

3.  INVERSE  SCATTERING  SOLUTION  STARTING  FROM  THE 
BODY  ITSELF.  (THIS  CASE  IS  A  SIMULATION  PROCEDURE, 
WHERE  THE  PROGRAM  RETRIEVES  BY  INVERSE  SCATTERING 
THE  SHAPE  OF  THE  KNOWN  BODY.  ( I NC*2  AND  INVER*l). 

THE  INPUTS  REQUIRED  ARE: 

SAME  AS  A  THRU  £  LISTED  ABOVE,PLUS  THE  LENGTH  OF 
A  SEGMENT  FOR  NEW  SUBDIVISION  OF  THE  BODY. (DSC) 

4.  INVERSE  SCATTERING  SOLUTION  STARTING  FROM  THE 
MEASURED  GIVEN  RAMP  RESPNSE  OF  AN  UNKNOWN  BODY  (HSFM). 
THIS  CASE  IS  THE  TRUE  COMPUTATION  OF  THE  SHAPE  OF  THE 
BODY  WHEN  ONLY  ITS  RAMP  RESPONSE  IS  KNOWN. 

(INC=2  AND  INVER»2) 

THE  INPUTS  REQUIRED  ARE: 

A.  THE  TIME  SAMPLES  OF  THE  RAMP  RESPONSE. (ENTERED 
IN  OATA  CAROS) 

B.  THE  LENGTH  OF  THE  DESIRED  SEGMENTS.  (OSO) 

C.  8  THRU  0  LISTED  ABOVE. 


NOTE:  NO  INVERSE  SCATTERING  IS  COMPUTED  WHEN  INC* I  . 
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parameter  definition: 

}f^^^*m******  ********** 


nC:  DEFINES  THE  TYPE  OF  INCIDENT  FIELD 

INC=l  THE  INCIDENT  FIELD  IS  A  GAUSSIAN  IMPULSE. 
INC=2  THE  INCIDENT  FIELD  IS  A  RAMP. 


IiNVER  :  DEFINES  THE  TYPE  OF  PROGRAM  REQUIRED 
INYER=0  NO  INVERSE  SCATTERING 

INVER=1  THE  PROGRAM  COMPUTES  THE  INVERSE  SCATTERING, 
STARTING  FROM  A  GIVEN  BODY  SHAPE. 

INVER=2  THE  PROGRAM  COMPUTES  THE  INVERSE  SCATTERING 
FROM  A  GIVEN  RAMP  RESPONSE  OF  THE  BODY. 


NS:  NUMBER  OF  SEGMENTS  SUBDIVIDING  THE  C3NTCUR  f<=2D) 

ITS:  NUMBER  OF  TIME  SAMPLES  NEEDED  I <=  60) 

or:  length  of  THE  TIME  INCREMENT  (SECONDS) 

0S(20)  :  length  OF  THE  SEGMENTS  SUBDIVIDING  THE 

CONTOUR  OF  THE  BODY. ( MAXI  MUM  20  SEGMENTS) 

TO  BE  ENTERO  IN  DATA  STATEMENT. 

Z(21):  2-CQORDINATE  OF  THE  SPACE  SAMPLE  POINTS  ON  THE 
CONTOUR  OF  THE  BODY. ( MAX  I  MUM  21  POINTS  ENTERED 
DELIMITING  20  SEGMENTS ).( METER) 

RJ  (21)  :RA0IUS  OF  THE  BODY  AT  THE  SPACE  SAMPLE  POINTS 
DEFINED  BY  Z.(MAX.  21  ). (METER) 

CJS(20,60):  SPACE  TIME  SAMPLES  OF  THE  S  COMPONENT  OF 
THE  CURRENTS  INDUCED  ON  THE  BODY.  (UP  TO 
20  SPACE  SAMPLES  AND  60  TIME  SAMPLES  EACH) 

CJPHI(20,60) ;  S’ACE  TIME  SAMPLES  OF  THE  PHI  COMPONENT 
OF  THE  CURRENTS  INDUCED  ON  THE  BODY,  (UP  TO 
20  SPACE  SAMPLES  AND  60  TIME  SAMPLES  EACH) 

HSF  (60):  TIME  SAMPLES  OF  THE  SCATTERED  FIELD. UP  TO 
60  SAMPLES. 

HSFM(60):  TIME  SAMPLES  OF  THE  INITIAL  RAMP  RESPONSE. 

(  THE  MEASURED  RAMP  RESPONSE) 

HSFC(60|:  TIME  SAMPLES  OF  THE  COMPUTED  RAMP  RESPONSE 

OF  THE  BODIES  OBTAINED  BY  INVERSE  SCATTERING 

HP (20):  NUMBER  OF  PATCHES  DN  THE  RINGS  DEFINED  BY 
EACH  SEGMENT. 

B(20);  ANGLE  SUSTAINED  BY  THE  PATCHES  ON  EACH  RING. 

ALFAK.;  SLOPE  OF  THE  CONTOUR  DF  THE  BODY  CONTOUR 
FUNCTION  , MEASURED  AT  THE  CENTER  OF  THE 
EACH  SEGMENT. 

SI  <20):  SIN(ALFAK)  FOR  EACH  SEGMENT. 

CO  (20):  COS( ALFAK)  FOR  EACH  SEGMENT. 

EPS:  NORMALIZED  MEAN  SQUARED  ERROR  ON  THE  BODY  SHAPE 
COMPUTED  BY  INVERSE  SCATTERING  PROCEDURE. 
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ITER:  NUMBER  OF  ITERATIONS  EXECUTED  FOR  INVERSE 
SCATTERING  SOLUTION. 

C:  SPEED  CF  LIGHT  (METER/SEC) 

C?  :  CONSTANT  C*DT. 

OTHER  PARAMETERS  USED  ARE  DEFINED  IN  THE  PROPER 
SUBROUTINE. 


REMARKS : 

1.  THE  SHAPE  OF  THE  BODY  IS  ENTERD  IN  DATA 
CAROS  WITH  THE  FORMAT  2F10. 5 

2.  THE  TIME  SAMPLES  OF  THE  MEASURED  RAMP 
RESPONSE  HSFM  ARE  ENTERD  IN  DATA  CAROS  WITH  THE 
FORMAT  1F10.5 

THE  VALUES  OF  THE  SAMPLES  MUST  BE  NORMAL IZEOt  THAT 
MEANS  MULTIPLIED  BY  C  AND  R  IRANGE  TO  TARGET) 

3.  IN  THIS  PROGRAM  THE  INCIDENT  FIELD  H  IS 
ASSUMED  TO  BE  VERTICAL  (HORIZONTAL  POLARIZATION). 
THE  DIFFERENCE  BETWEEN  THE  RESPONSE  TO  THE  TWO 
POLARIZATIONS  IS  ONLY  THE  SIGN  OF  THE  RESPONSE. 

(  SEE  ALSO  REMARKS  IN  SUBROUTINE  SHAPE) 


SUBROUTINES  USED: 

**^*$im*  ********** 

1.  SUBROUTINE  ‘FIELD*  :  COMPUTES  THE  BACK  SCATTERED 
FIELD  FOR  BOTH  RAMP  AND  IMPULSE  EXCITATION. 

2.  SUBROUTINE  'SHAPE*  :  COMPUTES  WITH  A  METHOD  OF 
ITERATIONS.  THE  SHAPE  CF  THE  BODY  FROM  ITS  RAMP 
RESPONSE. (SEE  SUBROUTINE  SHAPE  FOR  DETAILS.) 


3.  LIBRARY  SUBROUTINE  'PLOTP*  FOR  PLOTS  OF  FIELDS 
AND  SHAPE  OF  THE  BODY. 


METHOD  CF  SOLUTION: 

SEE  SUBROUTINE  FIELD  FOR  THE  COMPUTATION  OF  IMPULSE 
AND  RAMP  RESPONSES. 

SEE  SUBROUTINE  SHAPE  FOR  THE  COMPUTATION  OF  THE 
SHAPE  OF  THE  BODY. 

INVERSE  SCATTERING:  THE  FIRST  APPROXIMATION  OF  THE 
SHAPE  IS  OBTAINED  FROM  ITS  MEASURED  RAMP  RESPONSE. 

A  NEW  RAMP  RESPONSE  IS  COMPUTED  FOR  THAT  APPROXIMATE 
BODY. THE  PHYSICAL  OPTICS  PART  OF  THE  RESPONSE  IS 
EXTRACT  ED, AND  IS  USED  TO  COMPUTE  A  NEW  SHAPE, AND  SO 
ON  .UNTILL  A  MINIMUM  MEAN  SQUARED  ERROR  IS  OBTAINED 
BETWEEN  THE  LAST  COMPUTED  RESPONSE  AND  THE  MEASURED. 
(SEE  ALSO  SECTION  5  OF  THIS  THESIS) 
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c  *  * 

C  -  ^1AIN  PROGRAM  » 

C  *  * 

c 

C  SPHERE  EXAMPLE 

DIMENSION  X(6C) 

DIMENSION  CJS(20,60  I  , CJPH I ( 20 ,60 )  ,  riSF(60) 
DIMENSION  Z(21) ,R0< 21J ,NP(20) 

DIMENSION  SI  (20)  ,00  (20), 6  (20)  ,DS  (2D) 

DIMENSION  HSFM(60)  ,HSFC  (60) 

COMMON  HSF,Z,RO,ITS ,NS,NS1 

COMMON  XMAX 

Data  hsfm/60*d.o/ 

DATA  HSFC/60*0.D/ 

DATA  0S/20*0. 20900/ 

C 

C 

C  START  PROGRAM 

C 

C  DEFINE  TYPE  OF  PROGRAM, INVER  AND  INC 

r 

INC=2 

INVER=1 

C 

C  CQMPJTATION  PARAMETERS 

C 

NS=15 

NS1=NS>1 

ITS=60 

1T=ITS 

0SO=0.209 

CsS.OE+S 

0T=0.2/C 

C?=C*0T 

IF  (  INC.EQ. 1)  GO  TO  ll 
IF  ( INVER  .Eg. 2)  GO  TO  13 

r 

C  STARTING  POINT  FOR  INVER=1  AND  INVER=0 

C 

11  CONTINUE 

r 

C  READ  THE  SHAPE  OF  THE  BODY 

DO  10  R=l,NSi 
REA0(5,9510)  Z(K),RO(K» 

9513  FQRMAT(2F10.5) 

10  CONTINUE 

C 

C  PLOT  THE  INITIAL  SHAPE  OF  THE  BODY 

NS2=NSl-»-l 

XMAX  =  Z (NSL )*1.25 
Z{NS2)=  XMAX 
RQ  (NS2)=  Z(NS2I 
HRITE(6,9620) 

9620  FORMAT  (IHl) 

MQ0CUR=0 

CALL  PLOTP(  Z,  R0,NS2  ,;M00CUR) 
ftRlTE  (6,9621) 

9621  FORMAT  ( / / / , 20X , •  I  NIT  I AL  SHAPE  OF  BODY' , // t21X , 
2  'X-AXIS  Z,  Y-AXIS  RADIUS  (RO)') 

C 

c 
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c 

c 

c 


c 

c 

c 

c 


15 

c 


c 

c 

c 


12 


C 

19 

c 

c 

c 


910 

920 

C 

c 


c 

c 

c 

c 

c 


c 

c 

c 

c 


940 

c 

950 

C 

c 

c 

c 

c 


c 

c 

c 

c 


COMPUTE  THE  BACK- SCATTERED  FIELD 
CALL  FIELD  I0S,C?  ,  OT,  lUC  » 


END  OF  PROGRAM  FOR  INVER=0  AND  INVER*1 


IFi  INC.EQ.IJ  GO  TO  99 
IF  (INVER  .Eq.O)  GO  TO  99 
IF  (  INVER  .EO.l)  GO  TO  19 
CONTINUE 


STARTING  POINT  FOR  rNVER=2 
iNC  =  2 

READ  THE  GIVEN  RAMP  RESPONSE 


9520 


DO  12  N=l,ITS 

READ  (5,9520)  HSFM  (N) 

FORMAT  (IFIO.5) 

HSF (N)=HSFM(N) 

CONTINUE 
GO  TO  920 


CONTINUE 

STORE  THE  COMPUTED  RAMP  RESPONSE 


DO  910  N»1,ITS 
HSFM(NI*HSF<N) 


CONTINUE 

CONTINUE 


STARTING  POINT  FOR  INVERSE  SCATTERING  COMPUTATION 
ITER=1 

COMPUTATION  OF  THE  FIRST  AP PROX I MATI ON  OF  THE  SHAP 


CALL  SHAPE  ( OS , I TER,CP  .DSO , I T ) 


COMPUTE  SUM  OF  SQUARES  OF  THE  TIME  SAMPLES  OF  THE 
ORIGINAL  RAMP  RESPONSE. 


ITS=IT+10 
SIGMAl^O. 0 
JO  940  N=  1,  ITS 
SICMA1  =  SIGMA1«-HSFM(NI**2 
CONTINUE 


CONTINUE 

COMPUTE  THE  RESPONSE  OF  THE  APPROXIMATE  BODY. 


CALL  FIELD  ( 0  S,  CP  ,  DT,  INC  I 
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r 

c 


960 

C 

c 

c 

9666 


961 

C 

C 

C 


c 

c 

c 

c 

c 

c 

9649 


9651 

970 

C 

c 

c 

c 

c 

r 

w 

c 

c 

999 

c 

c 

c 

c 

9625 

99 

C 

r 

C 


CO'IPUTE  NOR.'IALIZED  MEAN  SOJARED  ERROR. 

SIG.MA2=0. 0 
03  960  N=lf ITS 

SIGMA2=SI  GMA2«-(HSFM(NI-HSF(N)  )**2 
CONTINUE 

EPS=SIGMA2/SIGMA1 
WRITE  ERROR 


WRITE  (6,9666)  EPS 

FORMAT  (  // ,15X, 'NORMALIZED  MEAN  SQ. 

IF  (ITER.GT.l)  GO  TO  961 
£3S1=1. 

CONTINUE 


ERROR*  ,1F10. 5) 


END  THE  PROGRAM  WHEN  THE  ERROR  IS  MINIMUM 

IF< EPS.GT.EPSl)  GO  TO  999 
ePSl=«EPS 

CONTINUE  WHEN  THE  ERROR  IS  NOT  YET  MINIMIZED 
ITER*ITER*-1 

EXTRACTION  OF  THE  PHYSICAL  OPTICS  RESPONSE 
WRITE  (6,9649) 

FORMAT  (IHl,  lOX, 'PHYSICAL  OPTICS  RESPONSE ',//, 9X ,* N* , 
2  10X,*H( N) • ,//) 

DO  970  N=1,ITS 
HSFC(N)*HSFC(N)+HSF(N) 

HSF(N)»ITER*HSFM(N)-HSFC(N) 

WRITE  (6,9651)  N,HSF(N) 

FQRMAT(5X,I5,5X,E12.5) 

CONTINUE 

USE  PHYSICAL  OPTICS  RESPONSE  TO  COMPUTE  A  NEW  SHAPE 
CALL  SHAPE  ( OS , I TER,CP  ,OSa, I T ) 

RETURN  TO  COMPUTE  THE  RESPONSE  TO  THE  NEW  30DY 

GO  TO  950 
CONTINUE 

END  OF  PROGRAM  WHEN  MINIMUM  ERROR  ACQUIRED 
WRITE  FINAL  ERROR 

WRITE  (6,9625)  ITER,EPS1 

FORMAT  (//,20X, 'FINAL  SHAPE  A FTER  • , 14, • I TERA T IONS • , /// 
2  ,20X. 'NORMALIZED  MEAN  SQUARED  ERROR  IS',1F10.5) 

CONTINUE 

END  OF  PROGRAM 

STOP 

END 


123 


ooooooonoooooooooooooooooooonooonnonoonooonooooorioo 


*4141^4141  «4i4i4i  41 «  4t4i4i4t4i 

*  4> 

♦  SUBROUTINE  FIELD  ♦ 

»  * 

«4>4i«4i««4i#  4:*  *41 4>4‘4>«  *  *4i<>#*4> 


PURPOSE  : 

«*«4t4i*4« 


THIS  SUBROUTINE  PROVIDES  THE  SOLUTION  TO  THE 
MAGNETIC  FIELD  INTEGRAL  EQUATION  APPLIED  TO  A 
BODY  OF  REVQLUTIONfEXCITEO  BY  A  GAUSSIAN  PULSE 
OR  A  RAMP  FIELDS. THE  CURRENTS  INDUCED  ON  THE  BODY 
AND  THE  BACK-SCATTERED  FIELD  ARE  COMPUTED. 

USAGE: 

«*»«4<« 


CALL  FIELD  IDS  #  CP  ,  DT,  INCI 


PARAMETER  DEFINITION: 

4i4^»*4<  4t*4(4i4i 


HI :  INCIDENT  FIELD 
C3,C4:  CONSTANTS. 


METHOD  OF  SOLUTION: 

4i4«)4  44‘4>  4‘4i4  4i4<4e4>**4  4i4< 


THE  METHOD  OF  SOLUTION  IS  DESCRIBED  IN  SECTION  4 
OF  THIS  THESIS. 


124 


or»o^a  r>oo  r.oo  ivjoooo  or>o  c'^oo  ooo  o  ono  non 


C 

c 


SJBROUTINE  FIELD  (  OS,CP  ,0T,  I\C  ) 


OIMEMSION  X(60J 

DIMENSION  CJS(20tt.O  )  ,  CJ  PH  IC  20  1 60  )  ,  HSF(60) 
DIMENSION  Z(21),R0(21) ,NP(20) 

DIMENSION  SI  (20»  t CO  ( 2C > , B  (20)  ,DS  (20) 

COMMON  HSF.Z.RO, ITSfNStNSl 

CONSTANT  DEFINITION 

PI=3.  141562654 
C=3 .oe+8 

C-»=-l./  (0T*4.  0*:) 

Ci=i./(  PI*2.0) 

GAUSSIAN  IMPULSE  PARAMETERS 

AA*  0  <6  E ♦ 9 
TMAX*0T*5.0 
A2  =  AA*AA 

DO  20  K=l,NS 
K1=K+1 

COMPUTE  the  slope  OF  THE  SEGMENTS 

WK*(R0(KI)-R0(K))/(Z(K1)-Z(K}  ) 

ALFAK=ATAN( WK J 

SI( K)=SIN(ALFAK) 

CC(  K)=CQS<ALFAK) 

compute  the  central  space  sample  points 

Z(K)*(Z(Kl)+Z(K) )/2.0 

R0(K)*(R0(K1)«>R0(K)»/2.0 

compute  the  NJM3ER  OF  PATCHES  IN  EACH  RING 

PR*  Pl*R0(K)/0S(K) 

NP( K) =2*1 FIX(  PK) 

IF  (NP(K) .GT.O)  GO  TO  23 
NP( R)=l 

COMPUTE  THE  ANGLE  SUSTAINED  BY  THE  PATCHES  OF  EACH 
R  ING. 

CONTINUE 

a(K)  =  2.0*Pl/NP(K) 

initial  CONDITIONS  OF  THE  CURRENTS 

TZK=1 .0+Z (K )/CP 
DO  25  N=l,2 

IF  (N.LT.TZK)  GO  TO  26 
IF(  I.NC.EQ.1  )  GO  TO  24 

WHEN  INCIDENT  FIELD  IS  A  RAMP 


HI=(N-1)  *CP-Z(i<) 
GO  TO  27 
continue 


WHEN  INCIDENT  FIELD  IS  AN  IMPULSE 


arg=(n-:) •0T'Tmax-zik)/c 

HI*(N-1)*£XP( -A2*ARG*ARGI 
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FCUR  POINT  LAGRANGI4N  INTERPOLATION 


Yi=T2*T3*T4*C JS( I . NQ2) /6.0-T3*T4*Tl*r JS( I , NQP )/2.0 
2  ♦T4*Tl*T2*CJS<  1 1  NQ I  /2. 0- Tl  ♦T2*  T3*CJS(  I  ,NO'<I)/6.0 

Y2=T2»T3*T4*C  JPHl  (  I  tNQ2  I/6.0-T3*T4*Tl*CJPHn  I , NQP i/ 2.0 
2  ♦T4*Tl*T2*CJPHI < I»NQ  1/2. 0-Tl*T2*T3*CJPHr ( I ,NQM»/6.0 
0JS=(T2*T3+r2*T4*T3*T4»*CJSI I,NQ2 1/6.0 

2  -( T3*T4+T4*TI+T1*T3I*CJSC I»NQP)/2.0 

3  ♦•(T4*Tl^T4*T2«-Tl*T2)*CJSt  l»Na»  /2.0 

4  -(Tl*T3+T2*T3+T2*Tl)*CJS< I,NQM)/6.0 
DJPHI=(T2^T3+T2*T4+T3*T4»*CJPHI ( I ,NC2 1/6.0 

2  -(T3*T4+T4<‘T1+T  l*T3)*CJPHn  I  , NQP  1/2.0 

3  ♦(T4*TI*T4*T2  +  T1*T2»*CJPHU  l,NQ)/2.0 

4  -(Tl*T3♦T2♦T3«•T2♦Tl)♦CJP^I^  ,NQM»/6.0 
GO  TO  440 

CONTINUE 

ThQ  point  LAGRANGIAN  INTERPOLATION 

Yl»T4*CjSU,Na  )-T3*CJS(IfNQMi 
Y2*T4*CJPHI ( I ,NQ  )-T3*CJPHI (I ,NQM) 

OJS=CJS(I#NQ  )-CJS(I,NQMI 
□JPHI»CJPHI  ( I  fNQ  )-CJPHlU#NQMI 
GO  TO  440 
CONTINUE 

THREE  POINT  LAGRANGIAN  INTERPOLATION 


Y1*T4*T3*CJS( I,NQP»/2.0-T2*T4*CJS(IfNQ  ) 

2  ♦T2*T3*CJS( I,NQMI/2.0  .  ^ 

r2  =  T4*T3*CJPHI(  I,NQPI/2.0-T2*T4#C  JPHKl  »NQ  ) 

2  ♦T2*T3*CJPHI( I.NQMI/2.0 
0JS  =  (T4+T3»  ♦CJS(  I  .NQP»/2.0-(T24-T4I*CJS(  ItNO) 

2  ♦(r2*T3)*CJS(  IfNOMI/2.0 

0JPHI=<T4-*-T3)  ♦CJPHI  (  I , NQP >/2. 0- I T 2>T4) *0 JPHI  (  I,NQ  I 
2  ♦|T2*T3)»CJPHI(I,NOMI/2.0 

GO  TO  440 
CONTINUE 

IF  IN.EQ.NQPI  go  TO  400 
IF  (N.EQ.NQ2)  GO  TO  433 

THREE  POINT  LAGRANGIAN  INTERPOLATION 

Y1*T2*T3«‘CJS(  I,NQ2l/2.0-T3*Tl*CJS(  I,NQP» 

2  ♦T2*Tl*CJS( IiNQ) /2.0 
Y2=«T2*T3*CJPHn  I,NQ2l/2.0-T3*Tl*CJPHI(  I,NQP» 

2  ♦T2*ti*CJPHI( I,NQI/2.0 
0JS  =  <T2*T3I*CJSU  iNQ2)/2.0-<T3*Tl  l♦CJS(  I  fNQP) 

2  ♦<T2+Tl »*CJS( I,NQI/2.0 

0JPHI»(T2+T3»  *CJPHn  I,NQ2l/2.0-(T3+Tl)*CJPHI(  I.NQPJ 
2  ♦(T2  +  TI  )*CJPHI{I,NQI/2.0 

GO  TO  440 
CONTINUE 

140  POINT  LAGRANGIAN  INTERPOLATION 

Yl»T3*CJS(  I  fNQPI-T2*CJSnfNQ) 

Y2»T3*CJPHI ( I ,NQP l-T2*CJPHI  ( I  ,NQI 
OJS*CJS(I fNQP »-CJS< I,NQ» 

0JPHI»CJPHI (ItNQPl-CJPHII If NQ) 
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c 

c 

c 


400 

C 

300 


C 

E 

C 

c 

c 

c 

c 


c 

c 

c 


240 

c 


250 

C 

C 

c 


s 

200 

C 

C 

c 

100 

c 

c 

c 


FINAL  COMPUTATION  OF  THE  INTEGRAL 
V=2l I i-li K) 

Fl*V^CM*SI(  n*(RQ(K)-R3(I  )*CN»*Ca(  I  ) 

Xl*Fl*CM 

X2*V*( 

U  =  V*SI(  I)*SI(KI*R0(KI*CQ4K)*S1(  I |-R0( !» *S I ( <1 ♦C0( I ) 
F2»(R0Ul  *CM-aO<  I )  )*CO|<)«-CM*V*SnK) 

Yll*Yl*Xl 
Y22»Y2*X2 
XPl*F2*CM 
XP2»U*C  SM*SM ) 

YPl*Y2*XPl 
YP2  =  Y  l*XP2 

Cl=OS(  n  *R0(  n*B(  I  l/tCP*R2» 

R3»R*R2 

C2*0S(  n*  ROI  I  ^♦a(I)/R3 
CJPHl  ((L,N)»CJPHUlL,NH>Cl*(DJPHI*XPl+0JS*XP2) 

2  +02*1  YPl«-YP2» 

CJSIK.nI*  CJS<K,N)*C1*(0JPHI*X2  ♦OJS*Xl  ) *02* ( Y1 !♦ Y22 ) 
CONTINUE 

CONTINUE 

C JS(K,NI»C JS(K»NI*C3 
CJPHK  K*NI  =CJPHI  (K,NI  ♦C3 


END  OF  COMPUTATION  OF  THE  INTEGRAL 


AOOITION  OF  THE  CURRENT  INDUCED  BY  THE  INCIDENT 
FIELD  AT  THE  COMPUTATION  POINT, (PHYSICAL  OPTICS  PART 

TZKal.0>2 (K)/CP 
T0*T2K«-2.*TMAX/0T 
IF  (N.LT.TZKI  GO  TO  200 
IF(  INC. £3.1)  GO  TO  240 

WHEN  INCIDENT  FIELD  IS  A  RAMP 

HI=IN-1) ♦CP-Z(<) 

GO  TO  250 
CONTINUE 

IF  (N.GT.TDI  GO  TO  200 

WHEN  INCIDENT  FIELD  IS  AN  IMPULSE 

ARG*(N-l)  ♦OT-TMAX-Z(KI/C 
HlaEXP(-A24ARG*ARG) 

CONTINUE 

FINAL  CURRENT  AT  THE  K-TH  SPACE  SAMPLE  AT  TIME  N 

CJS(K,NI«CJS(Kf NI-2.0*HI 

CJPHU  K,N|aCJPHI(K,N)-2.0*HI*SIlK) 

NEXT  SPACE  SAMPLE 

CONTINUE 

NEXT  TIME 

CONTINUE 

END  OF  COMPUTATION  OF  ALL  SPACE-TIME  CURRENTS 
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I 


5623 


110 

9723 

C 

c 

c 

c 

c 

111 

c 

c 

G 


C 

C 

c 

c 

c 

c 


c 

c 

c 

c 


c 

c 

c 

c 


602 

C 

8 

C 

603 


IF  (INC. 60.  1)  G3  TO  110 
•(RITE (6, 9623  » 

FORMAT! 131, lOX, 'NORMALIZED  FIELO  (H*C*R»',//, 

2  lOX,  'RAMP  RESPONSE* ,//,9X, 'N* , 10X,'H( N) • ,//) 

GO  TO  111 

CONTINUE 

NRITE(6,9723) 

FOfiMATI  HI  ,10X, 'NORMALIZED  FIELD  (H*R»  ',//, 

2  lOX,  'IMPULSE  RE SPONSE* ,//,9X,* N* ,10X,* H(N) • ,//l 


START  COMPUTATION  OF  THE  BACK-SCATTERED  FIELO 


CONTINUE 

AT  THE  N-TH  TIME 

00  500  N=l,  ITS 
HSF(NI*0.0 

IF  (N.EQ. 1)  GO  TO  650 
FROM  EACH  SAMPLE 
DO  600  K=1,NS 

COMPUTATION  OF  THE  RETARDATION  TIME 

a*N-z(K)/cp 

NQ=IFIX(QJ 

INTERPOLATION  OF  THE  TIME  DERIVATIVE  OF  THE  CURRENT 
AT  THE  RETARDED  TIME 


NaP*NQ*- 1 

NQMsNQ-  I 
N02*NQ+2 
T1=Q-NQ2 
T2=C-NQP 
T3=Q-NQ 
T4aQ-NQM 

IF  (Q.LE.l.OI  GO  TO  600 
IF  (Q.LE.2.01  GO  TO  601 

IF  (N.EQ.NQPJ  GO  TO  602 
IF  (N.EQ.NQ2I  GO  TO  603 


FOUR  POINT  LAGRANGIAN  INTERPOLATION 

DJS»<T2*T3+T2*T4+T3*T4I*CJS(K,NQ2 1/6.0 

2  -( T3*T4+T 4*Tl ♦Tl*T3l*CJ SI K,NQP 1/2.0 

3  ♦(T4*Tl+T4*T2+Tl*T2)*CJS(K,NQ)/2.0 

4  -(Tl«'T3*T2*T3*-T2*Tl)*CJSI  K,NQMI/6.0 
0JPHIa(T2*T3+T2*T4+T3*T4)*CJPHI(K,NQ2»/6.0 

2  -<T3*T4+T4*Tl+TI*T3»*CJPMl(K,NQPI/2.0 

3  ♦(T4*TUT4*T2  +  Tl*T2)*CJPHI(K,NQ»/2.0 

4  -I Tl*T3*T2*T3*T2*Tl»*CJPHI(K,NQM)/6.0 
GO  TO  62  0 

CONTINUE 


TWO  POINT  LAGRANGIAN  INTERPOLATION 

DJS  «  CJSiKfNQ  I-  CJSIK.NQMI 
OUPHI»CJPHI(K,NQ  »-CJPHlU,NQMI 
GO  TO  620 
CONTINUE 


I 

I 

I 
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THRSe  POINT  LAGRANGIAN  INTERPOLATION 

0JS  =  IT4+T3}  ♦CJS(K,NQP)/2.0-CT2<-T4)*CJS(K,NQ) 

Z  *IT2*T3»*t:jS(K,NQMI/2.0 

0JPHI=(T4+T3»  ♦CJPHI  (K,NQP)/2.0-(T2*-T4»»CjPHl  (  K,NQ  ) 
2  ♦(T2<-T3J*CJPHI(K»NQMI/2.0 

GO  TO  620 
CaNTiNOe 

I'r  (N.EQ.NaPJ  GC  TO  600 
IF  (N.ta.NQZ)  GO  TO  604 

THREE  POINT  LAGRANGIAN  INTERPOLATION 

0JS*(T2  +  T3)  *CJS(K,NQ2)/2.0-(T3*-Tl  )*CJS(K,NQP) 

2  ♦(T2«-Tl»*CJS(K,NQl/2.0 

0JPHI*(T2*T3I KaPHl (K,NQ2 I /2. 0- I T3+T1 I *0 JPHI ( K,NQP) 
2  ♦(T2^Tl  )*CJPHnK,NQ)/2.0 

GO  TO  62  0 
CONTINUE 

TWO  POINT  LAGRANGIAN  INTERPOLATION 

DJS  *  CJSIKtNQPI-  CJSfK»NQ» 

DJPHI»CJPHI (K,NQP  »-CJPHIlKtNQ) 

CONTINUE 

CONTRIBUTION  OF  THE  K-TH  SPACE  SAMPLE 

HSF(N)»HSF(  N)*ROIKl*(SI(KI*DJS  +  OJPHn*OSlK) 

NEXT  SPACE  SAMPLE 

CONTINUE 

CONTINUE 

FINAL  VALUE  OF  THE  N-TH  TIME  SAMPLE  OF  THE  FIELD 


9610 

C 


HSF(NI»C4*HSF(N  » 
X(N)=N*CP 

WRITE(6,9610I  NfHSFINI 
FORMAT! 5X, 15, 5X,E12.5» 

NEXT  TI<ME  SAMPLE 

CONTINUE 


9630 


9624 


9724 


END  OF  COMPUTATION  OF  THE  FIELD 

PLOT  BACK-SCATTEREO  FIELD  VERSUS  TIME! L IGHT-METERS) 

WRITE!6, 96301 
FORMAT(lHl) 

MODCUR*0 

CALL  PLOTPCX.HSF, ITS,MOOCUR» 

IF  (INC.EQ.lJ  GO  TO  112 
MRITE  16,96241 

FORMAT!// ♦14X.*  NQRMALUEO  BACK  SCATTERED  FIELD*,//, 

2  14X,*RAHP  ReSPQNSE',//jl4X, 

3  *X-AXIS  TIME  1C*0TP  ,10X,»Y-AXIS  FIELD  IH*C*R»'» 
GO  TO  113 

CONTINUE 
WRITE  16,9724) 

FORMAT!//, 14X,*  NORMALIZED  BACK  SCATTERED  FIELD',//, 

2  14X, 'IMPULSE  RESPONSE*, //,14X, 

3  U-AXIS  TIME  (C*On',lOX, 'Y-AXIS  FIELD  (H*RI*) 
CONTINUE 

RETURN 

END 
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Ai***1f*m**************^** 

«  * 

SJBRJUTIME  SHAPE  * 

*  * 

««««*«  ««  «**4i « 

PURPOSE; 

*«*«**«4i 

THIS  SUeROUTINE  COMPUTES  THE  SHAPE  CF  THE 
BODY  FROM  ITS  APPROXIMATE  COMPUTED  PHYSICAL-OPTICS 
RESPONSE. 

IT  FINDS  THE  CONTOUR  FUNCTION  RO(Z)  AND  DIVIDES  IT 
INTO  APPROXIMATELY  EQUAL  SEGMENTS  (OS). 

USAGE: 

**«*«* 

CALL  SHAPE  ( CS , I  TER , CP ,DSO ,  IT  ) 

THE  DATA  NEEDED  AT  THE  INPUT  ARE  OSOtANO  MSF. 

PARAMETER  DEFINITION: 

********************* 


ZMAX;  the  MAXIMUM  LENGTH  OF  THE  BODY  ON  THE  Z-AXIS. 
(CORRESPONDING  TO  THE  POINT  WHERE  THE  RAMP  RESPONSE 
CHANGES  ITS  SIGN) 

IT:  THE  TRANSIT  TIME  NEEDED  TO  OBTAIN  THE  REFLECTION 
FROM  THE  POINT  LOCATED  AT  ZMAX. (IT  IS  THE  NUMBER 
OF  TIME  INCREMENTS). 


DSO:  OPTIMAL  LENGTH  OF  A  SEGMENT ( NEEDED  FOR  THE 
COMPUTED  SHAPE  SUBDIVISION)  .(METER) 


ALL  OTHER  PARAMETERS  ALREADY  DEFINED  IN  MAIN  PROGRAM 


REMARKS  : 

******** 


IN  THIS  PROGRAM  THE  RAMP  RESPONSE  USED  ASSUMES  A 
VERTICAL  H  FIELD  INCIDENT  (HORIZONTAL  POLARIZATION) 
THIS  GIVES  A  POSITIVE  INITIAL  RESPONSE. 

FDR  VERTICAL  POLARIZATION  THE  INITIAL  RESPONSE  IS 
NEGATIVE.  (THE  WHOLE  RESPONSE  CHANGES  ITS  SIGN) 

SINCE  THE  RADIUS  OF  THE  BODY  AT  EACH  POINT  IS 
PROPORTIONAL  TO  SQRT  OF  THE  INITIAL  RAMP  RESPONSE, 

IT  MUST  BE  ENSURED  THAT  THE  RAMP  USED  HAS  PROPER  SIGN. 


METHOD  CF  SOLUTION: 


first  STEP  IS  TO  FIND  ZMAX, WHICH 
CORRESPONDS  TO  THE  POINT  WHERE  THE  RESPONSE 
CHANGES  ITS  SIGN. 

THE  CONTOUR  FUNCTION  OF  THE  BODY  RO(Z)  IS 
PROPOKriONAL  TO  THE  SQRT  OF  THE  PHYSICAL-OPTICS 
RAMP  RESPONSE. 

THE  BODY  IS  DIVIDED  INTO  APPROXIMATELY  EQUAL 
SEGMENTS, AND  THE  SPACE  SAMPLES  DELIMITING  THE 
SEGMENTS  ARE  FOUND  USING  A  'CUT  AND  TRY  METHOD* 
WHICH  CONVERGES  TO  THE  DESIRED  LENGTH  OF  SEGMENT. 
THE  RADIUS  AT  THE  SPACE  SAMPLE  POINTS  ARE 
INTERPOLATED  FROM  THE  KNOWN  SAMPLES  OF  THE  RAMP 
RESPONSE. 

THE  DATA  OBTAINED  Z ,RO,OS, NS, I S  USED  FOR  THE  NEXT 
COMPUTATION  OF  THE  FIELD. 
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c 

c 

SJBKQUTIME  SHAPE  ( DS , ITER , CP ,0  SO 1 1 T ) 

OIMENSION  OS(20)  tRQ<2U,  Z  (  2 1 )  ,  HSF  (  60 ) 

COMMON  HSF,ZtRQ»lTS»N5,NSl 
COMMON  XMAX 
C 

C  START  PROGRAM 

C 

C 

C  FIND  LENGTH  OF  BODY  ( ZMAX  ) AND  ITS  CORRESPONDING 

C  TIME  OF  APPEARANCE  IN  THE  RAMP  RESPONSE  (IT) 

DSP=OSG 

DO  710  N=l,  ITS 
IF(HSF(NJ  .LT.0.0)  GO  TO  720 
710  CONTINUE 

720  CONTINUE 

IT  =N-1 

AO=l./( l.-HSF(N)/HSF( IT  )l 
Z  UX=(  IT  -2>A0)*CP/2. 

C 

c 

C  COMPUTE  SHAPE  OF  BODY  UP  TO  ZMAX. 

C 

54  CONTINUE 
Z( 1  )=0.0 

RO  (1)=S3RT(2*HSF(2)  ) 

00  50  K=lflR 

K1=K+1 

'J  =  0SP/2.0 

55  CONTINUE 
NN*1 

C 

C  ADJUST  Z  TO  OBTAIN  PROPER  OS 

C 

Z(K1)=Z(K)+U 
53  CONTINUE 

IF(Z(K1).GT.ZMAXIG0  TO  56 
GO  TO  57 

56  CONTINUE 
Z(K1) =ZMAX 

RO(K1)=0.0 
GO  TO  51 

57  CONTINUE 

C  COMPUTE  CORRESPONDING  RADIUS  ( RO )  OF  THE  BODY 

C 

G=2*Z(K1)  /CP+2. 

M=IFIX(G) 

Ml*M+l 

R02  =  (HSF(  MH-HSF(M)  )♦(  G-Ml^-HSFl  M  ) 
R0(K1)=SQRT(R02*2.) 

c 

c  COMPUTE  LENGTH  OF  SEGMENT  OS 

C 

51  CONTINUE 

OS  (K)*SQRT(( Z(K1 )-Z(K) )**2+< R0( M  )-R0( X)  )**2) 
A=OS(K)/DSP 

S  ALL  OBTAINED  SEGMENTS  MUST  BE  WITHIN  1/100  OF  OSD 

C 

IF  (A.GT.l.Ol)  SO  TO  58 
1F( Z(K1 ). EQ.ZMAXIGO  TO  60 
IF  (A.LT.0.99)  GO  TO  58 
GO  TO  50 
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c 

C  ADJUST  z  AND  GO  BACK  TO  CHECK  FOR  OS 

C 

53  CONTINUE 

NN=NN*1 

IF  ( A.GT. l.OJ  GO  TO  90 
Z(K1»=Z(K1)+U/NN 
GJ  TO  53 
90  CONTINUE 

Z(K1)=Z<K1)-U/NN 
IF  (Z(K1 ) .GT.Z(K) I  GO  TO  53 
Z(K1)=Z(K1 ) VU/NN 
NN=NN+L 
GO  TO  90 
50  CONTINUE 

C 

C  NEXT  SEGMENT 

C 

60  CONTINUE 
NS1»K1 

^lS=NSl-l 

C 

C  ADJUST  FOR  LAST  SEGMENT 

C 

UZaOSINSJ  /OSP 
IF(U2.LT.O. 10 J  30  TO  62 
IF(U2,LT.0.95J  GO  TO  61 
GO  TO  70 

61  CONTINUE 
U3«(Z(NSll-Z( NSI )/( NS-l » 

NSlaNSl-l 

NS«NSl-i 
00  80  K»l,NS 
M«K*1 

Z(K1)=Z(  Kll  ♦K*'J3 
Z(NSII»ZMAX 
G^2«Z<K1> /CP+2. 

M*IFIX(G) 

R02*(HSF(Ml »-HSF(M»  !♦( 3-M )+HSF( M) 

RO(KUaS0Rr(RO2*2. ) 

RO(NS1)*0,0 

OS  (KI=SQRT(( Z(Kl »-Z(K»  » ♦♦2  +  ( R0( K 1 )-R0< K ) }**Z J 
33  CONTINUE 
GO  TO  70 

62  CONTINUE 
Z(NS)*ZMAX 
R0(NS)=0.0 
NS2=NS-l 

0S(NS2)aSQRT( (ZMAX-Z(NS2) »**2+R0(NS2)**2) 

NS1»NS1-1 
NS  *NSl-l 
73  CONTINUE 

C 

C  SHAPE  OF  BODY  COMPLETED  AND  SUBDIVIDED 

C  WRITE  NUMBER  OF  SEGMENTS  OBTAINED 

C 

WRITE  (6f 9601  )  ITER«NS 

9601  F0RMAT(lHl.5X,'NEW  NUMBER  OF  SEGMENTS  AFTER  ITERATION* 

2  fI4,'  IS*  ,l‘t,/7/.lOX,'Z*,lOX,*RO' ,10X,»0S' f/////> 

OS  (NS1)»99.9 
DO  40  N«LtNSl 

WRITE(6,9602)  Z( N) f ROI N) ,0S  (NI 

9602  FORMAT(5X»3F10.5I 

43  CONTINUE 

C 
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ooo 


Plot  the  cn.Nnu;^  function  of  the  body 

NS2=NS1>1 
Z(NS2)=  XMAX 
RO  (NS2»=  ZIMS2) 

WRITE(6,9603» 

9603  FORMAT!  IHU 
MQDCUR=0 

CALL  PL0TP( Zt R0|NS2fM00CUR) 

WRITE  (6,9604)  ITER 

9604  FORMAT!// ,20X, 'SHAPE  OF  BODY  AFTER  ITERATION  NUMBER 
2  I4,//,20X, 'X-AXIS  Z  Y-AXIS  RADIUS  !R0)») 

<1 

C  END  OF  subroutine  SHAPE 


riooooooo 


EXAMPLE  OF  OATA  CARDS 


0.0 

0.0 

0.  02183 

0.  20791 

0.08645 

0.  4067‘4 

0. 19098 

0. 58779 

0.33087 

0. 74314 

0.  50 

0. 86603 

0.69098 

0. 95106 

0. 89547 

0.  99452 

1.10453 

0.99452 

1. 30902 

0.  95106 

1.50 

0. 86603 

1.66913 

0.74314 

L. 80902 

0.  58779 

1.91355 

0. 40574 

1.97815 

0.20791 

2.0 

0.  0 
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APPENDIX  F 


RELATIONSHIP  BETWEEN  RAMP  RESPONSE  AND  SHAPE 


The  validity  of  (5.2)  will  be  demonstrated  analytically 
here  in  the  case  of  a  metallic  body  of  revolution.  We  first 
express  the  axially  directed  incident  ramp  field  as: 


H^(2,t)  = 


/  t  -  2/C  t  ^  z/C 

0  t  <  z/C 


(F.  1) 


The  physical  optics  current  induced  on  the  body  is  given  by: 


J  (r ,  t) 


2a  xHi(?,t) 
n  K 


(F.2) 


The  backscattered  far-zone  field  is  given  by  equation  (3.8) 
and  its  scalar  equivalent  for  the  body  of  revolution  is 
expressed  by  (3.16).  Substitution  of  the  current  in  (3.16) 
gives  the  physical  optics  ramp  response  as  follows: 


3 

“r^^^  =  “  4^^  i  |^[-4H^(z' ,T)sin  i-lp'ds*  (F.3) 

o  0 

The  time  derivative  of  the  incident  field  expressed  by  (F.l) 
is  unity: 


3t  at  ' 

Substitution  of  (F.4)  into  (F.3) 


t  <  z/C 

gives: 


(F.4) 
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H^(t) 


0 


J  sin  a*  c ' ds ' 


(F.5; 


From  the  geometry  of  the  body  of  revolution  we  note  that: 


do'  =  sin  a*  ds' 


Substitution  in  (F.5)  gives  the  following  expression: 


H^(t) 


r  C 
o  0 


P  (2) 

<  p ' dp • 


1  2  ,  , 
P  (z) 


2r  C 
o 


(F.6) 


Since  the  area  function  is  given  by: 

A (2)  =  r  (z) 

the  physical  optics  ramp  response  (F.6)  will  be  related  to 
A (z)  by 


H^(t) 


2nr  C 
o 


A(z) 


(F.7) 
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